! IMPORTANT NOTE:
! if not running within the GFDL FMS framework, comment out the following
! line:

#define GFDL_COMPATIBLE_MICROP

#ifdef GFDL_COMPATIBLE_MICROP
module micro_mg_mod
! this is ncar routine micro_mg1_5_CESM1.2.F90.
#else
module micro_mg1_5
#endif

!---------------------------------------------------------------------------------
! Purpose:
!   MG microphysics version 1.5 - Update of MG microphysics and jumping-off
!       point for the development of MG2
!
! Author: Andrew Gettelman, Hugh Morrison.
! Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan
! Version 2: Development begun: September 2011
! invoked in CAM by specifying -microphys=mg1.5
! 
! for questions contact Hugh Morrison, Andrew Gettelman
! e-mail: morrison@ucar.edu, andrew@ucar.edu
!
! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
! microphysics in cooperation with the MG liquid microphysics. This is
! controlled by the do_cldice variable.
!
! NOTE: If do_cldice is false, then MG microphysics should not update CLDICE
! or NUMICE; however, it is assumed that the other microphysics scheme will have
! updated CLDICE and NUMICE. The other microphysics should handle the following
! processes that would have been done by MG:
!   - Detrainment (liquid and ice)
!   - Homogeneous ice nucleation
!   - Heterogeneous ice nucleation
!   - Bergeron process
!   - Melting of ice
!   - Freezing of cloud drops
!   - Autoconversion (ice -> snow)
!   - Growth/Sublimation of ice
!   - Sedimentation of ice
!---------------------------------------------------------------------------------
! Based on micro_mg (restructuring of former cldwat2m_micro)
! Author: Andrew Gettelman, Hugh Morrison.
! Contributions from: Xiaohong Liu and Steve Ghan
! December 2005-May 2010
! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
!                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
! for questions contact Hugh Morrison, Andrew Gettelman
! e-mail: morrison@ucar.edu, andrew@ucar.edu
!---------------------------------------------------------------------------------

#ifdef GFDL_COMPATIBLE_MICROP
!--------------------------------------------------------------------------
! Interfaces, diagnostics, used modules, constants modified for use in GFDL
! based models
! Huan Guo and Rick Hemler, 2010 - 2013
!--------------------------------------------------------------------------
#endif

! Code comments added by HM, 093011
! General code structure:
!
! Code is divided into two main subroutines:
!   subroutine micro_mg_init --> initializes microphysics routine, should be called
!                                  once at start of simulation
!   subroutine micro_mg_tend --> main microphysics routine to be called each time step
!
! List of external functions:
!   qsat_water --> for calculating saturation vapor pressure with respect to liquid water
!   qsat_ice --> for calculating saturation vapor pressure with respect to ice
!   gamma   --> standard mathematical gamma function
! .........................................................................
! List of inputs through use statement in fortran90:
! Variable Name                      Description                Units
! .........................................................................
! gravit          acceleration due to gravity                    m s-2
! rair            dry air gas constant for air                  J kg-1 K-1
! tmelt           temperature of melting point for water          K
! cpair           specific heat at constant pressure for dry air J kg-1 K-1
! rh2o            gas constant for water vapor                  J kg-1 K-1
! latvap          latent heat of vaporization                   J kg-1
! latice          latent heat of fusion                         J kg-1
! qsat_water      external function for calculating liquid water
!                 saturation vapor pressure/humidity              -
! qsat_ice        external function for calculating ice
!                 saturation vapor pressure/humidity              pa
! rhmini          relative humidity threshold parameter for
!                 nucleating ice                                  -
! .........................................................................
! NOTE: List of all inputs/outputs passed through the call/subroutine statement
!       for micro_mg_tend is given below at the start of subroutine micro_mg_tend.
!---------------------------------------------------------------------------------

! Procedures required:
! 1) An implementation of the gamma function (if not intrinsic).
! 2) saturation vapor pressure and specific humidity over water
! 3) svp over ice

#ifdef GFDL_COMPATIBLE_MICROP
use gamma_mg_mod,              only: gamma =>gamma_mg
use lscloud_types_mod,         only: diag_id_type, diag_pt_type
use mpp_mod,                   only: input_nml_file
use fms_mod,                   only: mpp_pe, file_exist, error_mesg,  &
                                     open_namelist_file, FATAL, &
                                     stdlog, write_version_number, &
                                     check_nml_error, close_file, &
                                     mpp_root_pe
use simple_pdf_mod,            only: simple_pdf
use sat_vapor_pres_mod,        only: lookup_es2, lookup_es3, compute_qs
use physics_radiation_exch_mod, only : exchange_control_type


#else

#ifndef HAVE_GAMMA_INTRINSICS
use shr_spfn_mod, only: gamma => shr_spfn_gamma_nonintrinsic
#endif

use wv_sat_methods, only: &
     qsat_water => wv_sat_qsat_water, &
     qsat_ice => wv_sat_qsat_ice
#endif

implicit none
private
save

public :: &
     micro_mg_init, &
     micro_mg_get_cols, &
     micro_mg_tend

#ifdef GFDL_COMPATIBLE_MICROP
!------------------------------------------------------------------------
!--version number--------------------------------------------------------

character(len=128) :: Version = '$Id$'
character(len=128) :: Tagname = '$Name$'
#endif
 
#ifndef GFDL_COMPATIBLE_MICROP
! switch for specification rather than prediction of droplet and crystal number
! note: number will be adjusted as needed to keep mean size within bounds,
! even when specified droplet or ice number is used

! ***note: Even if constant cloud ice number is set, ice number is allowed
! to evolve based on process rates. This is needed in order to calculate
! the change in mass due to ice nucleation. All other ice microphysical
! processes are consistent with the specified constant ice number if
! this switch is turned on.

! nccons = .true. to specify constant cloud droplet number
! nicons = .true. to specify constant cloud ice number

logical, parameter, public :: nccons = .false.
logical, parameter, public :: nicons = .false.
#else
! nccons and nicons provided by micro_mg_nml. 
#endif

!=========================================================
! Private module parameters
!=========================================================

integer, parameter :: r8 = selected_real_kind(12)   ! 8 byte real

real(r8), parameter :: pi = 3.14159265358979323846_r8 ! To 20 digits; more than enough
                                                          ! to reach the limit of double precision

real(r8), parameter :: omsm   = 0.99999_r8    ! number near unity for round-off issues

#ifndef GFDL_COMPATIBLE_MICROP
! parameters for specified ice and droplet number concentration
! note: these are local in-cloud values, not grid-mean
real(r8), parameter :: ncnst = 100.e6_r8    ! droplet num concentration when nccons=.true. (m-3)
real(r8), parameter :: ninst = 0.1e6_r8     ! ice num concentration when nicons=.true. (m-3)
#else
! ncnst and ninst provided by micro_mg_nml.   
#endif

! rhow used to be passed in in init, but then overwritten.
! For now, just setting it as a parameter here
#ifndef GFDL_COMPATIBLE_MICROP
real(r8), parameter :: rhosn = 250._r8  ! bulk density snow
real(r8), parameter :: rhoi = 500._r8   ! bulk density ice
#else
!  rhosn, rhoi are in nml used by GFDL code. see below.
!  rhosn was added in nml because Marc Salzmann used a value of 100 
!  with his model, while MG use a value of 250., as seen here.
#endif
real(r8), parameter :: rhow = 1000._r8  ! bulk density liquid

! fall speed parameters, V = aD^b (V is in m/s)
! droplets
real(r8), parameter :: ac = 3.e7_r8
real(r8), parameter :: bc = 2._r8
! snow
real(r8), parameter :: as = 11.72_r8
real(r8), parameter :: bs = 0.41_r8
! cloud ice
real(r8), parameter :: ai = 700._r8
real(r8), parameter :: bi = 1._r8
! rain
real(r8), parameter :: ar = 841.99667_r8
real(r8), parameter :: br = 0.8_r8

! particle mass-diameter relationship
! currently we assume spherical particles for cloud ice/snow
! m = cD^d
! exponent
real(r8), parameter :: dsph = 3._r8

! ventilation parameters
! for snow
real(r8), parameter :: f1s = 0.86_r8
real(r8), parameter :: f2s = 0.28_r8
! for rain
real(r8), parameter :: f1r = 0.78_r8
real(r8), parameter :: f2r = 0.308_r8

! collection efficiencies
! aggregation of cloud ice and snow
real(r8), parameter :: eii = 0.1_r8

! autoconversion size threshold for cloud ice to snow (m)
#ifdef GFDL_COMPATIBLE_MICROP
! for gfdl, dcs is obtained from Exch_ctrl%(physics_driver_nml). 
! Old NCAR value is 400e-6, but GFDL default is 200.e-6 as 
! used by Marc Salzmann.
#else
real(r8), parameter :: dcs = 250.e-6_r8
#endif

! smallest mixing ratio considered in microphysics
real(r8), parameter :: qsmall = 1.e-18_r8  

! alternate threshold used for some in-cloud mmr
real(r8), parameter :: icsmall = 1.e-8_r8

! immersion freezing parameters, bigg 1953
real(r8), parameter :: bimm = 100._r8
real(r8), parameter :: aimm = 0.66_r8

#ifdef GFDL_COMPATIBLE_MICROP
! GFDL has min_diam_ice as a nml variable, so mi0 is calculated in
! the initialization subroutine using that value.
real(r8)            :: mi0 
#else
! mass of new crystal due to aerosol freezing and growth (kg)
real(r8), parameter :: mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)
real(r8) :: rhmini      ! Minimum rh for ice cloud fraction > 0.
logical :: activate_all_ice_always
logical :: do_hallet_mossop         
logical :: microp_uniform
logical :: do_cldice
#endif

!Range of cloudsat reflectivities (dBz) for analytic simulator
real(r8), parameter :: csmin = -30._r8
real(r8), parameter :: csmax = 26._r8
real(r8), parameter :: mindbz = -99._r8
real(r8), parameter :: minrefl = 1.26e-10_r8    ! minrefl = 10._r8**(mindbz/10._r8)

!=========================================================
! Constants set in initialization
!=========================================================

! Set using arguments to micro_mg_init
real(r8) :: g           ! gravity
real(r8) :: r           ! dry air gas constant
real(r8) :: rv          ! water vapor gas constant
real(r8) :: cpp         ! specific heat of dry air
real(r8) :: tmelt       ! freezing point of water (K)

! latent heats of:
real(r8) :: xxlv        ! vaporization
real(r8) :: xlf         ! freezing
real(r8) :: xxls        ! sublimation

!real(r8) :: rhmini      ! Minimum rh for ice cloud fraction > 0.

! flags
!logical :: microp_uniform
!logical :: do_cldice
!logical :: activate_all_ice_always
!logical :: do_hallet_mossop         

real(r8) :: rhosu       ! typical 850mn air density

real(r8) :: icenuct     ! ice nucleation temperature: currently -5 degrees C

real(r8) :: snowmelt    ! what temp to melt all snow: currently 2 degrees C
real(r8) :: rainfrze    ! what temp to freeze all rain: currently -5 degrees C

! additional constants to help speed up code
real(r8) :: cons1
real(r8) :: cons4
real(r8) :: cons5
real(r8) :: cons7
real(r8) :: cons8
real(r8) :: cons11
real(r8) :: cons13
real(r8) :: cons14
real(r8) :: cons16
real(r8) :: cons17
real(r8) :: cons22
real(r8) :: cons23
real(r8) :: cons24
real(r8) :: cons25
real(r8) :: cons27
real(r8) :: cons28

#ifdef GFDL_COMPATIBLE_MICROP
 
!real       :: dcs = 400.e-6_r8    !autoconversion size threshold
!real       :: min_diam_ice = 10.e-6_r8
real       :: dcs                 !autoconversion size threshold
real       :: min_diam_ice              
logical    :: allow_all_cldtop_collection = .false.
logical    :: rho_factor_in_max_vt = .true.
real       :: max_rho_factor_in_vt = 1.0
real       :: lowest_temp_for_sublimation = 180._r8
real       :: rhosn = 100._r8
real       :: rhoi  = 500._r8

! switch for specification rather than prediction of droplet and crystal number
! note: number will be adjusted as needed to keep mean size within bounds,
! even when specified droplet or ice number is used

! ***note: Even if constant cloud ice number is set, ice number is allowed
! to evolve based on process rates. This is needed in order to calculate
! the change in mass due to ice nucleation. All other ice microphysical
! processes are consistent with the specified constant ice number if
! this switch is turned on.

! nccons = .true. to specify constant cloud droplet number
! nicons = .true. to specify constant cloud ice number

logical :: nccons = .false.
logical :: nicons = .false.
 
! parameters for specified ice and droplet number concentration
! note: these are local in-cloud values, not grid-mean
real(r8) :: ncnst = 100.e6_r8    ! droplet num concentration when 
                                 ! nccons=.true. (m-3)
real(r8) :: ninst = 0.1e6_r8     ! ice num concentration when 
                                 ! nicons=.true. (m-3)

!--> cjg: modifications incorporated from Huan's code
logical    :: allow_rain_num_evap = .false.
real       :: accretion_scale = 1.0_r8

!<--cjg

logical           :: use_qcvar_in_accretion = .false.
real(r8), private :: qcvar_min4accr         = 0.1
real(r8), private :: qcvar_max4accr         = 0.5
real(r8), private :: accretion_scale_max    = 2.0
real(r8), private :: accr_scale

! <---h1g, 2012-06-12

logical           :: liu_in = .false.
                             ! True = Liu et al 2007 Ice nucleation
                             ! False = cooper fixed ice nucleation (MG2008)

logical           :: use_Meyers = .false.
! Ni (/m3) = 1000 * exp( (12.96* [(esl-esi)/esi]) - 0.639 )
!  Figure 9.3 of Rogers and Yau (1998) shows the nearly linear
!       variation of [(esl-esi)/esi] from 0. at 273.16K to 0.5 at 
!       233.16K.  Analytically this is parameterized as (tfreeze-T)/80.
!
!  Ni (/m3) = 1000 * exp( 12.96* (tfreeze-T)/80 - 0.639 )

logical           :: CldFallNew  = .true.
logical           :: SedBugFix   = .true.
logical           :: ActNew      = .true.
real              :: tc_cooper   = -35.0
real              :: IceFallFac  = 1.0
real              :: SnowFallFac = 1.0
logical           :: use_K2013_auto = .false.
logical           :: use_K2013_accr = .false.
logical           :: use_const_accretion_scale = .false.
logical           :: include_contact_freeze_in_berg = .false.
real              :: Bergeron_fac = 1.0
logical           :: include_homo_freez_in_qn = .true.      ! h1g, 2014-07-22
logical           :: include_contact_freez_in_qn = .true.   ! h1g, 2014-07-23
logical           :: include_immersion_freez_in_qn = .true. ! h1g, 2014-07-23
real              :: r_cri  = 0.0
logical           :: use_RK_auto = .false.
logical           :: use_RK_accr = .false.

logical           :: liq_sedi_mass_only = .false.   ! h1g, 2016-03-07
logical           :: no_liq_sedi        = .false.   ! h1g, 2016-03-07


real              ::  rhmini=0.80     ! minimum rh for ice cld fraction > 0
logical           :: activate_all_ice_always = .false.
logical           :: do_hallet_mossop = .false.
logical           ::  microp_uniform = .false.
                                      ! .true. = configure uniform for 
                                      ! sub-columns 
                                      ! .false. = use w/o sub-columns 
                                      ! (default) 
logical           ::  do_cldice = .true.
                                      ! .true. = do all processes (default)
                                      ! .false. = skip all processes 
                                      ! affecting cloud ice

!  <DATA NAME="activate_all_ice_always" TYPE="logical" DEFAULT=".true.">
!   should all potentially available ice particles be activated under all
!   conditions ? (or only when dqa is increasing)
!  </DATA>
!  <DATA NAME="do_hallet_mossop" TYPE="logical" DEFAULT=".false.">
!   the hallet-mossop process should be included in the NCAR microphysics?
!  </DATA>

namelist / micro_mg_nml /   &
                 rhmini, activate_all_ice_always, do_hallet_mossop, &
                  microp_uniform, do_cldice, &
                 allow_all_cldtop_collection, &
                 max_rho_factor_in_vt, &
                 rho_factor_in_max_vt, lowest_temp_for_sublimation, &
!--> cjg: modifications incorporated from Huan's code
                 allow_rain_num_evap, accretion_scale, &
!<--cjg
                 rhosn, rhoi,         &  ! h1g
                 nccons, ncnst,                &  ! cjg
                 nicons, ninst,                &  ! cjg
                 liu_in,                       &  ! h1g
                 use_qcvar_in_accretion,       &  ! h1g
                 qcvar_min4accr,               &  ! h1g
                 qcvar_max4accr,               &  ! h1g
                 accretion_scale_max,          &  ! h1g
                 CldFallNew, SedBugFix, ActNew, use_Meyers, tc_cooper, & 
                 IceFallFac, SnowFallFac, use_K2013_auto, use_K2013_accr, &
                 use_const_accretion_scale,   &
                 include_contact_freeze_in_berg, Bergeron_fac, &   
                 include_homo_freez_in_qn, include_contact_freez_in_qn, & 
                 include_immersion_freez_in_qn, r_cri, use_RK_auto,   &
                 use_RK_accr, liq_sedi_mass_only, no_liq_sedi



real              :: q_cri



#endif

! Generic interface for packing routines
interface pack_array
   module procedure pack_array_1Dr8
   module procedure pack_array_2Dr8
   module procedure pack_array_3Dr8
end interface

interface unpack_array
   module procedure unpack_array_1Dr8
   module procedure unpack_array_1Dr8_arrayfill
   module procedure unpack_array_2Dr8
   module procedure unpack_array_2Dr8_arrayfill
end interface

logical :: do_ice_nucl_wpdf

#ifdef GFDL_COMPATIBLE_MICROP
logical :: clubb_active
#endif
!===============================================================================
contains
!===============================================================================

#ifdef GFDL_COMPATIBLE_MICROP
subroutine micro_mg_init( &
     kind, gravit, rair, rh2o, cpair,    &
     tmelt_in, latvap, latice,           &
     do_ice_nucl_wpdf_in, errstring, Exch_ctrl)
#else
subroutine micro_mg_init( &
     kind, gravit, rair, rh2o, cpair,    &
     tmelt_in, latvap, latice,           &
     rhmini_in, microp_uniform_in, do_cldice_in,  &
     errstring)
#endif
  
  !----------------------------------------------------------------------- 
  ! 
  ! Purpose: 
  ! initialize constants for MG microphysics
  ! 
  ! Author: Andrew Gettelman Dec 2005
  ! 
  !-----------------------------------------------------------------------
#ifdef GFDL_COMPATIBLE_MICROP
  type(exchange_control_type), intent(in) :: Exch_ctrl   
#endif
  
  integer,  intent(in)  :: kind         ! Kind used for reals
  real(r8), intent(in)  :: gravit
  real(r8), intent(in)  :: rair
  real(r8), intent(in)  :: rh2o
  real(r8), intent(in)  :: cpair
  real(r8), intent(in)  :: tmelt_in     ! Freezing point of water (K)
  real(r8), intent(in)  :: latvap
  real(r8), intent(in)  :: latice
#ifndef GFDL_COMPATIBLE_MICROP
 real(r8), intent(in)  :: rhmini_in    ! Minimum rh for ice cloud fraction > 0.

 logical,  intent(in)  :: microp_uniform_in    ! .true. = configure uniform for sub-columns 
                                            ! .false. = use w/o sub-columns (standard)
 logical,  intent(in)  :: do_cldice_in     ! .true. = do all processes (standard)
                                            ! .false. = skip all processes affecting
                                            !           cloud ice
#endif

! logical, intent(in) :: activate_all_ice_always_in
! logical, intent(in) :: do_hallet_mossop_in
#ifdef GFDL_COMPATIBLE_MICROP
  logical, intent(in) :: do_ice_nucl_wpdf_in
#endif
  character(128), intent(out) :: errstring    ! Output status (non-blank for error return)
 
#ifdef GFDL_COMPATIBLE_MICROP
   INTEGER   :: unit, io, ierr, logunit
#endif
  !-----------------------------------------------------------------------

  errstring = ' '

  if( kind .ne. r8 ) then
     errstring = 'micro_mg_init: KIND of reals does not match'
     return
  endif

  ! declarations for MG code (transforms variable names)

  g= gravit                 ! gravity
  r= rair                   ! dry air gas constant: note units(phys_constants are in J/K/kmol)
  rv= rh2o                  ! water vapor gas constant
  cpp = cpair               ! specific heat of dry air
  tmelt = tmelt_in
#ifndef GFDL_COMPATIBLE_MICROP
  rhmini = rhmini_in
  microp_uniform = microp_uniform_in
  do_cldice  = do_cldice_in
#endif

  ! latent heats

  xxlv = latvap         ! latent heat vaporization
  xlf  = latice         ! latent heat freezing
  xxls = xxlv + xlf     ! latent heat of sublimation

  ! flags
! microp_uniform = microp_uniform_in
! do_cldice  = do_cldice_in
! activate_all_ice_always = activate_all_ice_always_in
! do_hallet_mossop = do_hallet_mossop_in

#ifdef GFDL_COMPATIBLE_MICROP

!---------------------------------------------------------------
!     process namelist
!---------------------------------------------------------------
#ifdef INTERNAL_FILE_NML
      read (input_nml_file, nml=micro_mg_nml, iostat=io)
      ierr = check_nml_error(io,'micro_mg_nml')
#else
      if ( file_exist('input.nml')) then
        unit = open_namelist_file ()
        ierr=1; do while (ierr /= 0)
        read  (unit, nml=micro_mg_nml, iostat=io, end=10)
        ierr = check_nml_error(io,'micro_mg_nml')
        enddo
10      call close_file (unit)
      endif
#endif

!-----------------------------------------------------------------------
!    write version and namelist to stdlog.
!-----------------------------------------------------------------------
      call write_version_number (version, tagname)
      logunit = stdlog()
      if (mpp_pe() == mpp_root_pe()) &
                      write (logunit, nml=micro_mg_nml)

       dcs = Exch_ctrl%dcs
       min_diam_ice = Exch_ctrl%min_diam_ice
   do_ice_nucl_wpdf = do_ice_nucl_wpdf_in
      
#endif

#ifdef GFDL_COMPATIBLE_MICROP
   mi0 = 4._r8/3._r8*pi*rhoi*(min_diam_ice)*(min_diam_ice)*(min_diam_ice)
   clubb_active=(Exch_ctrl%do_clubb>0)
!-----------------------------------------------------------------------
!    check for nml variable compatibility.
!-----------------------------------------------------------------------
   if ( liq_sedi_mass_only .and. no_liq_sedi        ) &
        call error_mesg ('micro_mg_mod',  &
           'liq_sedi_mass_only and no_liq_sedi cannot be true at the same time', FATAL)

#endif
  ! typical air density at 850 mb

  rhosu = 85000._r8/(rair * tmelt)

  ! Maximum temperature at which snow is allowed to exist
  snowmelt = tmelt + 2._r8
  ! Minimum temperature at which rain is allowed to exist
  rainfrze = tmelt - 5._r8

  ! Ice nucleation temperature
  icenuct  = tmelt - 5._r8

  ! Define constants to help speed up code (this limits calls to gamma function)
  ! Unused names: cons6, cons15, cons21, cons26
  cons1=gamma(1._r8+dsph)
  cons4=gamma(1._r8+br)
  cons5=gamma(4._r8+br)
  cons7=gamma(1._r8+bs)     
  cons8=gamma(4._r8+bs)     
  cons11=gamma(3._r8+bs)
  cons13=gamma(5._r8/2._r8+br/2._r8)
  cons14=gamma(5._r8/2._r8+bs/2._r8)
  cons16=gamma(1._r8+bi)
  cons17=gamma(4._r8+bi)
  cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
  cons23=dcs**3
  cons24=dcs**2
  cons25=dcs**bs
  cons27=xxlv**2
  cons28=xxls**2

end subroutine micro_mg_init

!===============================================================================
!microphysics routine for each timestep goes here...

#ifdef GFDL_COMPATIBLE_MICROP
subroutine micro_mg_tend ( &
      dqa_activation, total_activation, tiedtke_macrophysics, &
      j, jdim, mgncol,   mgcols, nlev, top_lev,  deltatin,           &
     tn,  qn, qcn, qin, ncn, nin, relvarn, accre_enhann,          &     
     pn,                 pdeln,              pint,               &
     cldn,               liqcldf,            icecldf,            &
     delta_cfin, D_eros_l4in, nerosc4in, D_eros_i4in, nerosi4in,   &
     dqcdtin, dqidtin, naain,    npccnin,  rndstn,   naconin,  &
     tlato,    qvlato,   qctendo,  qitendo,  nctendo,  nitendo,  &
     precto,   precio,  qsouto, rflxo,    sflxo,    &
     qrouto,   reff_raino,         reff_snowo,         &
     tnd_qsnown,         tnd_nsnown,         re_icen,            &
     errstring,  f_snow_berg,                                       &
     ssat_disposalout, n_diag_4d, diag_4d, diag_id, diag_pt        )

  !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
  ! e-mail: morrison@ucar.edu, andrew@ucar.edu

  ! input arguments
  logical, intent (in) :: dqa_activation
  logical, intent (in) :: total_activation
  logical, intent (in) :: tiedtke_macrophysics
  integer,  intent(in) :: j, jdim
  integer,  intent(in) :: mgncol           ! number of microphysics columns
  integer,  intent(in) :: mgcols(:)        ! list of microphysics columns
  integer,  intent(in) :: nlev             ! number of layers
  integer,  intent(in) :: top_lev          ! top level to do microphysics
  real(r8), intent(in) :: deltatin         ! time step (s)
  real(r8), intent(in) :: tn(:,:)          ! input temperature (K)
  real(r8), intent(in) :: qn(:,:)          ! input h20 vapor mixing ratio (kg/kg)
  real(r8), intent(in) :: relvarn(:,:)     ! relative variance of cloud water (-)
  real(r8), intent(in) :: accre_enhann(:,:)! optional accretion enhancement factor (-)

  ! note: all input cloud variables are grid-averaged
  real(r8), intent(in) :: qcn(:,:)       ! cloud water mixing ratio (kg/kg)
  real(r8), intent(in) :: qin(:,:)       ! cloud ice mixing ratio (kg/kg)
  real(r8), intent(in) :: ncn(:,:)       ! cloud water number conc (1/kg)
  real(r8), intent(in) :: nin(:,:)       ! cloud ice number conc (1/kg)
  real(r8), intent(in) :: pn(:,:)         ! air pressure (pa)
  real(r8), intent(in) :: pdeln(:,:)      ! pressure difference across level (pa)
  ! hm add 11-16-11, interface pressure
  real(r8), intent(in) :: pint(:,:)    ! level interface pressure (pa)
  real(r8), intent(in) :: cldn(:,:)      ! cloud fraction (no units)
  real(r8), intent(in) :: liqcldf(:,:)   ! liquid cloud fraction (no units)
  real(r8), intent(in) :: icecldf(:,:)   ! ice cloud fraction (no units)
  real(r8), intent(in) :: delta_cfin(:,:)        ! 
  real(r8), intent(inout) :: D_eros_l4in(:,:)           
  real(r8), intent(in) :: nerosc4in (:,:)            
  real(r8), intent(inout) :: D_eros_i4in(:,:)            
  real(r8), intent(in) :: nerosi4in (:,:)            
  real(r8), intent(in) :: dqcdtin   (:,:)            
  real(r8), intent(in) :: dqidtin   (:,:)        
  ! used for scavenging
  ! Inputs for aerosol activation
  real(r8), intent(in) :: naain(:,:)     ! ice nucleation number (from microp_aero_ts) (1/kg)
  real(r8), intent(in) :: npccnin(:,:)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)

  ! Note that for these variables, the dust bin is assumed to be the last index.
  ! (For example, in CAM, the last dimension is always size 4.)
  real(r8), intent(in) :: rndstn(:,:,:)  ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m)
  real(r8), intent(in) :: naconin(:,:,:) ! number in each dust bin, for contact freezing  (from microp_aero_ts) (1/m^3)

  ! Used with CARMA cirrus microphysics
  ! (or similar external microphysics model)
  real(r8), intent(in) :: tnd_qsnown(:,:) ! snow mass tendency (kg/kg/s)
  real(r8), intent(in) :: tnd_nsnown(:,:) ! snow number tendency (#/kg/s)
  real(r8), intent(in) :: re_icen(:,:)    ! ice effective radius (m)

  ! output arguments

  real(r8), intent(out) :: tlato(:,:)   ! latent heating rate       (W/kg)
  real(r8), intent(out) :: qvlato(:,:)  ! microphysical tendency qv (1/s)
  real(r8), intent(out) :: qctendo(:,:) ! microphysical tendency qc (1/s)
  real(r8), intent(out) :: qitendo(:,:) ! microphysical tendency qi (1/s)
  real(r8), intent(out) :: nctendo(:,:) ! microphysical tendency nc (1/(kg*s))
  real(r8), intent(out) :: nitendo(:,:) ! microphysical tendency ni (1/(kg*s))
  real(r8), intent(out) :: precto(:)    ! surface precip rate (m/s)
  real(r8), intent(out) :: precio(:)    ! cloud ice/snow precip rate (m/s)
  real(r8), intent(out) :: qsouto(:,:)  ! snow mixing ratio (kg/kg)
  real(r8), intent(out) :: rflxo(:,:)   ! grid-box average rain flux (kg m^-2 s^-1)
  real(r8), intent(out) :: sflxo(:,:)   ! grid-box average snow flux (kg m^-2 s^-1)
  real(r8), intent(out) :: qrouto(:,:)  ! grid-box average rain mixing ratio (kg/kg)
  real(r8), intent(out) :: reff_raino(:,:) ! rain effective radius (micron)
  real(r8), intent(out) :: reff_snowo(:,:) ! snow effective radius (micron)

  character(128),   intent(out) :: errstring        ! output status (non-blank for error return)
  real(r8), intent(out) :: f_snow_berg  (:,:) ! ratio of bergeron 
                                                     ! production of qi to 
                                                     ! sum of bergeron, 
                                                     ! riming and freezing
  real(r8), intent(out) :: ssat_disposalout(:,:)
                                 ! disposition of supersaturation at end 
                                 ! of step; 0.= no ssat, 1.= liq, 2.=ice)
INTEGER,INTENT(IN) :: n_diag_4d
REAL, dimension( :,:,:, 0: ), INTENT(INOUT) ::  diag_4d
TYPE(diag_id_type),INTENT(IN) :: diag_id
TYPE(diag_pt_type),INTENT(INout) :: diag_pt

#else
subroutine micro_mg_tend ( &
     mgncol,   mgcols,   nlev,     top_lev,  deltatin,           &
     tn,                 qn,                                     &
     qcn,                          qin,                          &
     ncn,                          nin,                          &
     relvarn,            accre_enhann,                           &     
     pn,                 pdeln,              pint,               &
     cldn,               liqcldf,            icecldf,            &
     rate1ord_cw2pr_st,  naain,    npccnin,  rndstn,   naconin,  &
     tlato,    qvlato,   qctendo,  qitendo,  nctendo,  nitendo,  &
     effco,    effco_fn, effio,              precto,   precio,   &
     nevapro, evapsnowo, praino,  prodsnowo, cmeouto,  deffio,   &
     pgamrado, lamcrado, qsouto,   dsouto,   rflxo,    sflxo,    &
     qrouto,             reff_raino,         reff_snowo,         &
     qcsevapo, qisevapo, qvreso,   cmeiout,  vtrmco,   vtrmio,   &
     qcsedteno,qisedteno,prao,     prco,     mnuccco,  mnuccto,  &
     msacwio,  psacwso,  bergso,   bergo,    melto,    homoo,    &
     qcreso,             prcio,    praio,    qireso,             &
     mnuccro,  pracso,   meltsdto, frzrdto,  mnuccdo,            &
     nrouto,   nsouto,   reflo,    areflo,   areflzo,  freflo,   &
     csrflo,   acsrflo,  fcsrflo,            rercldo,            &
     ncaio,    ncalo,    qrouto2,  qsouto2,  nrouto2,  nsouto2,  &
     drouto2,  dsouto2,  freqso,   freqro,   nficeo,             &
     tnd_qsnown,         tnd_nsnown,         re_icen,            &
     errstring)

  !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
  ! e-mail: morrison@ucar.edu, andrew@ucar.edu

  ! input arguments
  integer,  intent(in) :: mgncol                ! number of microphysics columns
  integer,  intent(in) :: mgcols(:)             ! list of microphysics columns
  integer,  intent(in) :: nlev                  ! number of layers
  integer,  intent(in) :: top_lev               ! top level to do microphysics
  real(r8), intent(in) :: deltatin              ! time step (s)
  real(r8), intent(in) :: tn(:,:)               ! input temperature (K)
  real(r8), intent(in) :: qn(:,:)               ! input h20 vapor mixing ratio (kg/kg)
  real(r8), intent(in) :: relvarn(:,:)          ! relative variance of cloud water (-)
  real(r8), intent(in) :: accre_enhann(:,:)     ! optional accretion enhancement factor (-)

  ! note: all input cloud variables are grid-averaged
  real(r8), intent(in) :: qcn(:,:)       ! cloud water mixing ratio (kg/kg)
  real(r8), intent(in) :: qin(:,:)       ! cloud ice mixing ratio (kg/kg)
  real(r8), intent(in) :: ncn(:,:)       ! cloud water number conc (1/kg)
  real(r8), intent(in) :: nin(:,:)       ! cloud ice number conc (1/kg)
  real(r8), intent(in) :: pn(:,:)         ! air pressure (pa)
  real(r8), intent(in) :: pdeln(:,:)      ! pressure difference across level (pa)
  ! hm add 11-16-11, interface pressure
  real(r8), intent(in) :: pint(:,:)    ! level interface pressure (pa)
  real(r8), intent(in) :: cldn(:,:)      ! cloud fraction (no units)
  real(r8), intent(in) :: liqcldf(:,:)   ! liquid cloud fraction (no units)
  real(r8), intent(in) :: icecldf(:,:)   ! ice cloud fraction (no units)
  ! used for scavenging
  ! Inputs for aerosol activation
  real(r8), intent(in) :: naain(:,:)     ! ice nucleation number (from microp_aero_ts) (1/kg)
  real(r8), intent(in) :: npccnin(:,:)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)

  ! Note that for these variables, the dust bin is assumed to be the last index.
  ! (For example, in CAM, the last dimension is always size 4.)
  real(r8), intent(in) :: rndstn(:,:,:)  ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m)
  real(r8), intent(in) :: naconin(:,:,:) ! number in each dust bin, for contact freezing  (from microp_aero_ts) (1/m^3)

  ! Used with CARMA cirrus microphysics
  ! (or similar external microphysics model)
  real(r8), intent(in) :: tnd_qsnown(:,:) ! snow mass tendency (kg/kg/s)
  real(r8), intent(in) :: tnd_nsnown(:,:) ! snow number tendency (#/kg/s)
  real(r8), intent(in) :: re_icen(:,:)    ! ice effective radius (m)

  ! output arguments

  real(r8), intent(out) :: rate1ord_cw2pr_st(:,:)    ! 1st order rate for
  ! direct cw to precip conversion
  real(r8), intent(out) :: tlato(:,:)         ! latent heating rate       (W/kg)
  real(r8), intent(out) :: qvlato(:,:)        ! microphysical tendency qv (1/s)
  real(r8), intent(out) :: qctendo(:,:)       ! microphysical tendency qc (1/s)
  real(r8), intent(out) :: qitendo(:,:)       ! microphysical tendency qi (1/s)
  real(r8), intent(out) :: nctendo(:,:)       ! microphysical tendency nc (1/(kg*s))
  real(r8), intent(out) :: nitendo(:,:)       ! microphysical tendency ni (1/(kg*s))
  real(r8), intent(out) :: effco(:,:)         ! droplet effective radius (micron)
  real(r8), intent(out) :: effco_fn(:,:)      ! droplet effective radius, assuming nc = 1.e8 kg-1
  real(r8), intent(out) :: effio(:,:)         ! cloud ice effective radius (micron)
  real(r8), intent(out) :: precto(:)          ! surface precip rate (m/s)
  real(r8), intent(out) :: precio(:)          ! cloud ice/snow precip rate (m/s)
  real(r8), intent(out) :: nevapro(:,:)       ! evaporation rate of rain + snow (1/s)
  real(r8), intent(out) :: evapsnowo(:,:)     ! sublimation rate of snow (1/s)
  real(r8), intent(out) :: praino(:,:)        ! production of rain + snow (1/s)
  real(r8), intent(out) :: prodsnowo(:,:)     ! production of snow (1/s)
  real(r8), intent(out) :: cmeouto(:,:)       ! evap/sub of cloud (1/s)
  real(r8), intent(out) :: deffio(:,:)        ! ice effective diameter for optics (radiation) (micron)
  real(r8), intent(out) :: pgamrado(:,:)      ! ice gamma parameter for optics (radiation) (no units)
  real(r8), intent(out) :: lamcrado(:,:)      ! slope of droplet distribution for optics (radiation) (1/m)
  real(r8), intent(out) :: qsouto(:,:)        ! snow mixing ratio (kg/kg)
  real(r8), intent(out) :: dsouto(:,:)        ! snow diameter (m)
  real(r8), intent(out) :: rflxo(:,:)         ! grid-box average rain flux (kg m^-2 s^-1)
  real(r8), intent(out) :: sflxo(:,:)         ! grid-box average snow flux (kg m^-2 s^-1)
  real(r8), intent(out) :: qrouto(:,:)        ! grid-box average rain mixing ratio (kg/kg)
  real(r8), intent(out) :: reff_raino(:,:)    ! rain effective radius (micron)
  real(r8), intent(out) :: reff_snowo(:,:)    ! snow effective radius (micron)
  real(r8), intent(out) :: qcsevapo(:,:)      ! cloud water evaporation due to sedimentation (1/s)
  real(r8), intent(out) :: qisevapo(:,:)      ! cloud ice sublimation due to sublimation (1/s)
  real(r8), intent(out) :: qvreso(:,:)        ! residual condensation term to ensure RH < 100% (1/s)
  real(r8), intent(out) :: cmeiout(:,:)       ! grid-mean cloud ice sub/dep (1/s)
  real(r8), intent(out) :: vtrmco(:,:)        ! mass-weighted cloud water fallspeed (m/s)
  real(r8), intent(out) :: vtrmio(:,:)        ! mass-weighted cloud ice fallspeed (m/s)
  real(r8), intent(out) :: qcsedteno(:,:)     ! qc sedimentation tendency (1/s)
  real(r8), intent(out) :: qisedteno(:,:)     ! qi sedimentation tendency (1/s)

  ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
  real(r8), intent(out) :: prao(:,:)         ! accretion of cloud by rain 
  real(r8), intent(out) :: prco(:,:)         ! autoconversion of cloud to rain
  real(r8), intent(out) :: mnuccco(:,:)      ! mixing ratio tend due to immersion freezing
  real(r8), intent(out) :: mnuccto(:,:)      ! mixing ratio tend due to contact freezing
  real(r8), intent(out) :: msacwio(:,:)      ! mixing ratio tend due to H-M splintering
  real(r8), intent(out) :: psacwso(:,:)      ! collection of cloud water by snow
  real(r8), intent(out) :: bergso(:,:)       ! bergeron process on snow
  real(r8), intent(out) :: bergo(:,:)        ! bergeron process on cloud ice
  real(r8), intent(out) :: melto(:,:)        ! melting of cloud ice
  real(r8), intent(out) :: homoo(:,:)        ! homogeneous freezing cloud water
  real(r8), intent(out) :: qcreso(:,:)       ! residual cloud condensation due to removal of excess supersat
  real(r8), intent(out) :: prcio(:,:)        ! autoconversion of cloud ice to snow
  real(r8), intent(out) :: praio(:,:)        ! accretion of cloud ice by snow
  real(r8), intent(out) :: qireso(:,:)       ! residual ice deposition due to removal of excess supersat
  real(r8), intent(out) :: mnuccro(:,:)      ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
  real(r8), intent(out) :: pracso(:,:)       ! mixing ratio tendency due to accretion of rain by snow (1/s)
  real(r8), intent(out) :: meltsdto(:,:)     ! latent heating rate due to melting of snow  (W/kg)
  real(r8), intent(out) :: frzrdto(:,:)      ! latent heating rate due to homogeneous freezing of rain (W/kg)
  real(r8), intent(out) :: mnuccdo(:,:)      ! mass tendency from ice nucleation
  real(r8), intent(out) :: nrouto(:,:)        ! rain number concentration (1/m3)
  real(r8), intent(out) :: nsouto(:,:)        ! snow number concentration (1/m3)
  real(r8), intent(out) :: reflo(:,:)         ! analytic radar reflectivity        
  real(r8), intent(out) :: areflo(:,:)        ! average reflectivity will zero points outside valid range
  real(r8), intent(out) :: areflzo(:,:)       ! average reflectivity in z.
  real(r8), intent(out) :: freflo(:,:)        ! fractional occurrence of radar reflectivity
  real(r8), intent(out) :: csrflo(:,:)        ! cloudsat reflectivity 
  real(r8), intent(out) :: acsrflo(:,:)       ! cloudsat average
  real(r8), intent(out) :: fcsrflo(:,:)       ! cloudsat fractional occurrence of radar reflectivity
  real(r8), intent(out) :: rercldo(:,:)       ! effective radius calculation for rain + cloud
  real(r8), intent(out) :: ncaio(:,:)        ! output number conc of ice nuclei available (1/m3)
  real(r8), intent(out) :: ncalo(:,:)        ! output number conc of CCN (1/m3)
  real(r8), intent(out) :: qrouto2(:,:)       ! copy of qrout as used to compute drout2
  real(r8), intent(out) :: qsouto2(:,:)       ! copy of qsout as used to compute dsout2
  real(r8), intent(out) :: nrouto2(:,:)       ! copy of nrout as used to compute drout2
  real(r8), intent(out) :: nsouto2(:,:)       ! copy of nsout as used to compute dsout2
  real(r8), intent(out) :: drouto2(:,:)       ! mean rain particle diameter (m)
  real(r8), intent(out) :: dsouto2(:,:)       ! mean snow particle diameter (m)
  real(r8), intent(out) :: freqso(:,:)        ! fractional occurrence of snow
  real(r8), intent(out) :: freqro(:,:)        ! fractional occurrence of rain
  real(r8), intent(out) :: nficeo(:,:)        ! fractional occurrence of ice

  character(128),   intent(out) :: errstring        ! output status (non-blank for error return)
#endif

#ifdef GFDL_COMPATIBLE_MICROP
!   these variables are output by NCAR, but not GFDL, so need to be 
!   declared as local for GFDL implementation

  real(r8) :: diag_4l(mgncol,jdim,nlev,0:n_diag_4d) ! truncated diagnostic
                                                  ! array
  real(r8) :: f_snow_berg_l    (mgncol,nlev)   
  real(r8) :: delta_cf         (mgncol,nlev)   
  real(r8) :: D_eros_l4        (mgncol,nlev)   
  real(r8) :: nerosc4          (mgncol,nlev)   
  real(r8) :: D_eros_i4        (mgncol,nlev)   
  real(r8) :: nerosi4          (mgncol,nlev)   
  real(r8) :: dqcdt            (mgncol,nlev)   
  real(r8) :: dqidt            (mgncol,nlev)   
  real(r8) :: ssat_disposal    (mgncol,nlev)   
  real(r8) :: rate1ord_cw2pr_st(mgncol,nlev)    ! 1st order rate for
  ! direct cw to precip conversion
  real(r8) :: effco(mgncol,nlev)         ! droplet effective radius (micron)
  real(r8) :: effco_fn(mgncol,nlev)      ! droplet effective radius, assuming nc = 1.e8 kg-1
  real(r8) :: effio(mgncol, nlev)         ! cloud ice effective radius (micron)
  real(r8) :: nevapro(mgncol,nlev)       ! evaporation rate of rain + snow (1/s)
  real(r8) :: evapsnowo(mgncol,nlev)     ! sublimation rate of snow (1/s)
  real(r8) :: praino(mgncol,nlev)        ! production of rain + snow (1/s)
  real(r8) :: prodsnowo(mgncol,nlev)     ! production of snow (1/s)
  real(r8) :: cmeouto(mgncol,nlev)       ! evap/sub of cloud (1/s)
  real(r8) :: deffio(mgncol,nlev)        ! ice effective diameter for optics (radiation) (micron)
  real(r8) :: pgamrado(mgncol,nlev)      ! ice gamma parameter for optics (radiation) (no units)
  real(r8) :: lamcrado(mgncol,nlev)      ! slope of droplet distribution for optics (radiation) (1/m)
  real(r8) :: dsouto(mgncol,nlev)  ! snow diameter (m)
  real(r8) :: qcsevapo(mgncol,nlev) ! cloud water evaporation due to sedimentation (1/s)
  real(r8) :: qisevapo(mgncol,nlev)      ! cloud ice sublimation due to sublimation (1/s)
  real(r8) :: qvreso(mgncol,nlev)        ! residual condensation term to ensure RH < 100% (1/s)
  real(r8) :: cmeiout(mgncol,nlev)      ! grid-mean cloud ice sub/dep (1/s)
  real(r8) :: vtrmco(mgncol,nlev)        ! mass-weighted cloud water fallspeed (m/s)
  real(r8) :: vtrmio(mgncol,nlev)        ! mass-weighted cloud ice fallspeed (m/s)
  real(r8) :: qcsedteno(mgncol,nlev)     ! qc sedimentation tendency (1/s)
  real(r8) :: qisedteno(mgncol,nlev)     ! qi sedimentation tendency (1/s)

  ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
  real(r8) :: prao(mgncol,nlev)         ! accretion of cloud by rain 
  real(r8):: prco(mgncol,nlev)         ! autoconversion of cloud to rain
  real(r8) :: mnuccco(mgncol,nlev)      ! mixing ratio tend due to immersion freezing
  real(r8) :: mnuccto(mgncol,nlev)      ! mixing ratio tend due to contact freezing
  real(r8) :: msacwio(mgncol,nlev)      ! mixing ratio tend due to H-M splintering
  real(r8):: psacwso(mgncol,nlev)      ! collection of cloud water by snow
  real(r8) :: bergso(mgncol,nlev)       ! bergeron process on snow
  real(r8) :: bergo(mgncol,nlev)        ! bergeron process on cloud ice
  real(r8) :: melto(mgncol,nlev)        ! melting of cloud ice
  real(r8) :: homoo(mgncol,nlev)       ! homogeneous freezing cloud water
  real(r8) :: qcreso(mgncol,nlev)       ! residual cloud condensation due to removal of excess supersat
  real(r8) :: prcio(mgncol,nlev)        ! autoconversion of cloud ice to snow
  real(r8):: praio(mgncol,nlev)        ! accretion of cloud ice by snow
  real(r8) :: qireso(mgncol,nlev)       ! residual ice deposition due to removal of excess supersat
  real(r8) :: mnuccro(mgncol,nlev)      ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
  real(r8) :: pracso(mgncol,nlev)       ! mixing ratio tendency due to accretion of rain by snow (1/s)
  real(r8) :: meltsdto(mgncol,nlev)     ! latent heating rate due to melting of snow  (W/kg)
  real(r8) :: frzrdto(mgncol,nlev)      ! latent heating rate due to homogeneous freezing of rain (W/kg)
  real(r8) :: mnuccdo(mgncol,nlev)      ! mass tendency from ice nucleation
  real(r8) :: nrouto(mgncol,nlev)        ! rain number concentration (1/m3)
  real(r8) :: nsouto(mgncol,nlev)        ! snow number concentration (1/m3)
  real(r8) :: reflo(mgncol,nlev)         ! analytic radar reflectivity        
  real(r8) :: areflo(mgncol,nlev)        ! average reflectivity will zero points outside valid range
  real(r8) :: areflzo(mgncol,nlev)       ! average reflectivity in z.
  real(r8) :: freflo(mgncol,nlev)        ! fractional occurrence of radar reflectivity
  real(r8) :: csrflo(mgncol,nlev)        ! cloudsat reflectivity 
  real(r8) :: acsrflo(mgncol,nlev)       ! cloudsat average
  real(r8) :: fcsrflo(mgncol,nlev)       ! cloudsat fractional occurrence of radar reflectivity
  real(r8) :: rercldo(mgncol,nlev)       ! effective radius calculation for rain + cloud
  real(r8) :: ncaio(mgncol,nlev)        ! output number conc of ice nuclei available (1/m3)
  real(r8) :: ncalo(mgncol,nlev)        ! output number conc of CCN (1/m3)
  real(r8) :: qrouto2(mgncol,nlev)       ! copy of qrout as used to compute drout2
  real(r8) :: qsouto2(mgncol,nlev)       ! copy of qsout as used to compute dsout2
  real(r8) :: nrouto2(mgncol,nlev)       ! copy of nrout as used to compute drout2
  real(r8) :: nsouto2(mgncol,nlev)       ! copy of nsout as used to compute dsout2
  real(r8) :: drouto2(mgncol,nlev)       ! mean rain particle diameter (m)
  real(r8) :: dsouto2(mgncol,nlev)       ! mean snow particle diameter (m)
  real(r8) :: freqso(mgncol,nlev)        ! fractional occurrence of snow
  real(r8) :: freqro(mgncol,nlev)        ! fractional occurrence of rain
  real(r8):: nficeo(mgncol,nlev)        ! fractional occurrence of ice

!  these variables are only used in the GFDL implementation

  real(r8) :: cmelo(mgncol,nlev)   ! liquid condensation           
  real(r8) :: eroslo(mgncol,nlev)  ! liquid erosion                
  real(r8) :: erosio(mgncol,nlev)  ! ice erosion                
  real(r8) :: preo(mgncol,nlev)    ! rain evaporation 
  real(r8) :: prdso(mgncol,nlev)   ! snow sublimation 
  real(r8) :: npccn2(mgncol,nlev)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
  logical  :: do_berg1
  logical  :: limit_berg = .false.
  real(r8) :: berg_lim = 1.0e-6_r8
  real(r8) :: dum3                ! temporary dummy variable

! droplet number
   real(r8) :: nucclim(mgncol,nlev)
   real(r8) :: nucclimo(mgncol,nlev)
  real(r8) :: npccno(mgncol,nlev)
   real(r8) :: nnuccco(mgncol,nlev)
  real(r8) :: nnuccto(mgncol,nlev)
  real(r8) :: npsacwso(mgncol,nlev)
   real(r8) :: nsubco(mgncol,nlev)
   real(r8) :: nerosco(mgncol,nlev)
   real(r8) :: nprao(mgncol,nlev)
   real(r8) :: nprc1o(mgncol,nlev)
   real(r8) :: nerosc(mgncol,nlev)
   real(r8) :: D_eros_l(mgncol,nlev)

! cloud ice number
  real(r8) :: nucclim1i(mgncol,nlev)
  real(r8) :: nucclim1io(mgncol,nlev)
   real(r8) :: nnuccdo(mgncol,nlev)
   real(r8) :: nsacwio(mgncol,nlev)
 real(r8) :: nsubio(mgncol,nlev)
  real(r8) :: nerosio(mgncol,nlev)
  real(r8) :: nprcio(mgncol,nlev)
  real(r8) :: npraio(mgncol,nlev)
  real(r8) :: nerosi(mgncol,nlev)
 real(r8) :: D_eros_i(mgncol,nlev)


  real(r8) :: cmel     (mgncol,nlev)
  real(r8) :: cmel_orig(mgncol,nlev)
 real(r8) :: cmei_orig(mgncol,nlev)
  real(r8) :: berg_orig(mgncol,nlev)

 real(r8) :: sum_freeze(mgncol,nlev)
 real(r8) :: sum_freeze2(mgncol,nlev)
 real(r8) :: sum_rime  (mgncol,nlev)
 real(r8) :: sum_berg  (mgncol,nlev)
  real(r8) :: sum_ice_adj(mgncol,nlev)
  real(r8) :: sum_bergs (mgncol,nlev)
   real(r8) :: sum_cond  (mgncol,nlev)
   real(r8) :: sum_splinter(mgncol,nlev)
   real(r8) :: qldt_sum
    real(r8) :: eslt, esit, rhi, qs_d, tc
    real(r8) :: qs2d(mgncol,nlev)
   real(r8) :: qtot(mgncol,nlev)


#endif





  ! local workspace
  ! all units mks unless otherwise stated

  ! parameters
  real(r8), parameter :: mincld = 0.0001_r8     ! minimum allowed cloud fraction
  real(r8), parameter :: cdnl   = 0.e6_r8       ! cloud droplet number limiter
  ! assign number of sub-steps to iter
  ! use 2 sub-steps, following tests described in MG2008

  integer, parameter :: iter = 2

  ! local copies of input variables
  real(r8) :: q(mgncol,nlev)           ! water vapor mixing ratio (kg/kg)
  real(r8) :: t(mgncol,nlev)           ! temperature (K)
  real(r8) :: qc(mgncol,nlev)          ! cloud liquid mixing ratio (kg/kg)
  real(r8) :: qi(mgncol,nlev)          ! cloud ice mixing ratio (kg/kg)
  real(r8) :: nc(mgncol,nlev)          ! cloud liquid number concentration (1/kg)
  real(r8) :: ni(mgncol,nlev)          ! cloud liquid number concentration (1/kg)
  real(r8) :: p(mgncol,nlev)           ! pressure (Pa)
  real(r8) :: pdel(mgncol,nlev)        ! pressure difference across level (Pa)
  real(r8) :: relvar(mgncol,nlev)      ! relative variance of cloud water (-)
  real(r8) :: accre_enhan(mgncol,nlev) ! optional accretion enhancement factor (-)

  real(r8) :: naai(mgncol,nlev)    ! ice nucleation number (from microp_aero_ts) (1/kg)
  real(r8) :: npccn(mgncol,nlev)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)

  real(r8), allocatable :: rndst(:,:,:)
  real(r8), allocatable :: nacon(:,:,:)

  real(r8) :: tnd_qsnow(mgncol,nlev) ! snow mass tendency (kg/kg/s)
  real(r8) :: tnd_nsnow(mgncol,nlev) ! snow number tendency (#/kg/s)
  real(r8) :: re_ice(mgncol,nlev)    ! ice effective radius (m)

  ! Packed copies of output variables
  real(r8) :: tlat(mgncol,nlev)         ! latent heating rate       (W/kg)
  real(r8) :: qvlat(mgncol,nlev)        ! microphysical tendency qv (1/s)
  real(r8) :: qctend(mgncol,nlev)       ! microphysical tendency qc (1/s)
  real(r8) :: qitend(mgncol,nlev)       ! microphysical tendency qi (1/s)
  real(r8) :: nctend(mgncol,nlev)       ! microphysical tendency nc (1/(kg*s))
  real(r8) :: nitend(mgncol,nlev)       ! microphysical tendency ni (1/(kg*s))

  real(r8) :: effc(mgncol,nlev)         ! droplet effective radius (micron)
  real(r8) :: effc_fn(mgncol,nlev)      ! droplet effective radius, assuming nc = 1.e8 kg-1
  real(r8) :: effi(mgncol,nlev)         ! cloud ice effective radius (micron)

  real(r8) :: prect(mgncol)          ! surface precip rate (m/s)
  real(r8) :: preci(mgncol)          ! cloud ice/snow precip rate (m/s)

  real(r8) :: nevapr(mgncol,nlev)       ! evaporation rate of rain + snow (1/s)
  real(r8) :: evapsnow(mgncol,nlev)     ! sublimation rate of snow (1/s)
  real(r8) :: prain(mgncol,nlev)        ! production of rain + snow (1/s)
  real(r8) :: prodsnow(mgncol,nlev)     ! production of snow (1/s)
  real(r8) :: cmeout(mgncol,nlev)       ! evap/sub of cloud (1/s)
  real(r8) :: deffi(mgncol,nlev)        ! ice effective diameter for optics (radiation) (micron)
  real(r8) :: pgamrad(mgncol,nlev)      ! ice gamma parameter for optics (radiation) (no units)
  real(r8) :: lamcrad(mgncol,nlev)      ! slope of droplet distribution for optics (radiation) (1/m)


  real(r8) :: qsout(mgncol,nlev)        ! snow mixing ratio (kg/kg)
  real(r8) :: qsout2(mgncol,nlev)       ! copy of qsout as used to compute dsout2
  real(r8) :: nsout(mgncol,nlev)        ! snow number concentration (1/m3)
  real(r8) :: nsout2(mgncol,nlev)       ! copy of nsout as used to compute dsout2
  real(r8) :: dsout(mgncol,nlev)        ! snow diameter (m)
  real(r8) :: dsout2(mgncol,nlev)       ! mean snow particle diameter (m)

  real(r8) :: qrout(mgncol,nlev)        ! grid-box average rain mixing ratio (kg/kg)
  real(r8) :: qrout2(mgncol,nlev)       ! copy of qrout as used to compute drout2
  real(r8) :: nrout(mgncol,nlev)        ! rain number concentration (1/m3)
  real(r8) :: nrout2(mgncol,nlev)       ! copy of nrout as used to compute drout2
  real(r8) :: drout2(mgncol,nlev)       ! mean rain particle diameter (m)

  real(r8) :: reff_rain(mgncol,nlev)    ! rain effective radius (micron)
  real(r8) :: reff_snow(mgncol,nlev)    ! snow effective radius (micron)

  real(r8) :: freqs(mgncol,nlev)        ! fractional occurrence of snow
  real(r8) :: freqr(mgncol,nlev)        ! fractional occurrence of rain

  real(r8) :: rflx(mgncol,nlev+1)       ! grid-box average rain flux (kg m^-2 s^-1)
  real(r8) :: sflx(mgncol,nlev+1)       ! grid-box average snow flux (kg m^-2 s^-1)

  real(r8) :: qcsevap(mgncol,nlev)      ! cloud water evaporation due to sedimentation (1/s)
  real(r8) :: qisevap(mgncol,nlev)      ! cloud ice sublimation due to sublimation (1/s)
  real(r8) :: qvres(mgncol,nlev)        ! residual condensation term to ensure RH < 100% (1/s)
  real(r8) :: cmeitot(mgncol,nlev)      ! grid-mean cloud ice sub/dep (1/s)
  real(r8) :: vtrmc(mgncol,nlev)        ! mass-weighted cloud water fallspeed (m/s)
  real(r8) :: vtrmi(mgncol,nlev)        ! mass-weighted cloud ice fallspeed (m/s)
  real(r8) :: qcsedten(mgncol,nlev)     ! qc sedimentation tendency (1/s)
  real(r8) :: qisedten(mgncol,nlev)     ! qi sedimentation tendency (1/s)

  real(r8) :: pratot(mgncol,nlev)         ! accretion of cloud by rain 
  real(r8) :: prctot(mgncol,nlev)         ! autoconversion of cloud to rain
  real(r8) :: mnuccctot(mgncol,nlev)      ! mixing ratio tend due to immersion freezing
  real(r8) :: mnuccttot(mgncol,nlev)      ! mixing ratio tend due to contact freezing
  real(r8) :: msacwitot(mgncol,nlev)      ! mixing ratio tend due to H-M splintering
  real(r8) :: psacwstot(mgncol,nlev)      ! collection of cloud water by snow
  real(r8) :: bergstot(mgncol,nlev)       ! bergeron process on snow
  real(r8) :: bergtot(mgncol,nlev)        ! bergeron process on cloud ice
  real(r8) :: melttot(mgncol,nlev)        ! melting of cloud ice
  real(r8) :: homotot(mgncol,nlev)        ! homogeneous freezing cloud water
  real(r8) :: qcrestot(mgncol,nlev)       ! residual cloud condensation due to removal of excess supersat
  real(r8) :: prcitot(mgncol,nlev)        ! autoconversion of cloud ice to snow
  real(r8) :: praitot(mgncol,nlev)        ! accretion of cloud ice by snow
  real(r8) :: qirestot(mgncol,nlev)       ! residual ice deposition due to removal of excess supersat
  real(r8) :: mnuccrtot(mgncol,nlev)      ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
  real(r8) :: pracstot(mgncol,nlev)      ! mixing ratio tendency due to accretion of rain by snow (1/s)
  real(r8) :: mnuccdtot(mgncol,nlev)      ! mass tendency from ice nucleation
  real(r8) :: meltsdttot(mgncol,nlev)      ! latent heating rate due to melting of snow  (W/kg)
  real(r8) :: frzrdttot(mgncol,nlev)      ! latent heating rate due to homogeneous freezing of rain (W/kg)

  real(r8) :: refl(mgncol,nlev)         ! analytic radar reflectivity        
  real(r8) :: arefl(mgncol,nlev)        ! average reflectivity will zero points outside valid range
  real(r8) :: areflz(mgncol,nlev)       ! average reflectivity in z.
  real(r8) :: frefl(mgncol,nlev)        ! fractional occurrence of radar reflectivity
  real(r8) :: csrfl(mgncol,nlev)        ! cloudsat reflectivity 
  real(r8) :: acsrfl(mgncol,nlev)       ! cloudsat average
  real(r8) :: fcsrfl(mgncol,nlev)       ! cloudsat fractional occurrence of radar reflectivity

  real(r8) :: rercld(mgncol,nlev)       ! effective radius calculation for rain + cloud

  real(r8) :: nfice(mgncol,nlev)        ! fractional occurrence of ice

  real(r8) :: ncai(mgncol,nlev)   ! output number conc of ice nuclei available (1/m3)
  real(r8) :: ncal(mgncol,nlev)    ! output number conc of CCN (1/m3)

  ! general purpose variables
  real(r8) :: deltat            ! sub-time step (s)
  real(r8) :: mtime             ! the assumed ice nucleation timescale

  real(r8) :: dz(mgncol,nlev)         ! height difference across model vertical level
  real(r8) :: int_to_mid(mgncol, nlev) ! Coefficients for linear interpolation from
                                       ! interface to mid-level

  ! temporary variables for sub-stepping
  real(r8) :: tlat1(mgncol,nlev)
  real(r8) :: qvlat1(mgncol,nlev)
  real(r8) :: qctend1(mgncol,nlev)
  real(r8) :: qitend1(mgncol,nlev)
  real(r8) :: nctend1(mgncol,nlev)
  real(r8) :: nitend1(mgncol,nlev)
  real(r8) :: prect1(mgncol)
  real(r8) :: preci1(mgncol)

  ! physical properties of the air at a given point
  real(r8) :: rho(mgncol,nlev)    ! density (kg m-3)
  real(r8) :: dv(mgncol,nlev)     ! diffusivity of water vapor
  real(r8) :: mu(mgncol,nlev)     ! viscosity
  real(r8) :: sc(mgncol,nlev)     ! schmidt number
  real(r8) :: rhof(mgncol,nlev)   ! density correction factor for fallspeed

  ! cloud fractions
  real(r8) :: cldmax(mgncol,nlev) ! precip fraction assuming maximum overlap
  real(r8) :: cldm(mgncol,nlev)   ! cloud fraction
  real(r8) :: icldm(mgncol,nlev)  ! ice cloud fraction
  real(r8) :: lcldm(mgncol,nlev)  ! liq cloud fraction

  ! mass mixing ratios
  real(r8) :: qcic(mgncol,nlev)   ! in-cloud cloud liquid
  real(r8) :: qiic(mgncol,nlev)   ! in-cloud cloud ice
  real(r8) :: qsic(mgncol,nlev)   ! in-precip snow
  real(r8) :: qric(mgncol,nlev)   ! in-precip rain

  ! number concentrations
  real(r8) :: ncic(mgncol,nlev)   ! in-cloud droplet
  real(r8) :: niic(mgncol,nlev)   ! in-cloud cloud ice
  real(r8) :: nsic(mgncol,nlev)   ! in-precip snow
  real(r8) :: nric(mgncol,nlev)   ! in-precip rain
  ! maximum allowed ni value
  real(r8) :: nimax(mgncol,nlev)
#ifdef GFDL_COMPATIBLE_MICROP
  ! maximum allowed nc value
  real(r8) :: ncmax(mgncol,nlev)
#endif

  ! Size distribution parameters for:
  ! cloud ice
  real(r8) :: lami(mgncol,nlev)   ! slope
  real(r8) :: n0i(mgncol,nlev)    ! intercept
  ! cloud liquid
  real(r8) :: lamc(mgncol,nlev)   ! slope
  real(r8) :: pgam(mgncol,nlev)   ! spectral width parameter
  real(r8) :: cdist1(mgncol,nlev) ! droplet freezing calculation
  ! snow
  real(r8) :: lams(mgncol,nlev)   ! slope
  real(r8) :: n0s(mgncol,nlev)    ! intercept
  ! rain
  real(r8) :: lamr(mgncol,nlev)   ! slope
  real(r8) :: n0r(mgncol,nlev)    ! intercept

  ! combined size of precip & cloud drops
  integer :: arcld(mgncol, nlev)  ! averaging control flag

  ! Rates/tendencies due to:
  ! deposition/sublimation of cloud ice
  real(r8) :: cmei(mgncol,nlev)
  ! ice nucleation
  real(r8) :: nnuccd(mgncol,nlev) ! number rate from deposition/cond.-freezing
  real(r8) :: mnuccd(mgncol,nlev) ! mass mixing ratio
  ! freezing of cloud water
  real(r8) :: mnuccc(mgncol,nlev) ! mass mixing ratio
  real(r8) :: nnuccc(mgncol,nlev) ! number concentration
  ! contact freezing of cloud water
  real(r8) :: mnucct(mgncol,nlev) ! mass mixing ratio
  real(r8) :: nnucct(mgncol,nlev) ! number concentration
  ! HM ice multiplication
  real(r8) :: msacwi(mgncol,nlev) ! mass mixing ratio
  real(r8) :: nsacwi(mgncol,nlev) ! number conc
  ! autoconversion of cloud droplets
  real(r8) :: prc(mgncol,nlev)    ! mass mixing ratio
  real(r8) :: nprc(mgncol,nlev)   ! number concentration (rain)
  real(r8) :: nprc1(mgncol,nlev)  ! number concentration (cloud droplets)
  ! self-aggregation of snow
  real(r8) :: nsagg(mgncol,nlev)  ! number concentration
  ! self-collection of rain
  real(r8) :: nragg(mgncol,nlev)  ! number concentration
  ! collection of droplets by snow
  real(r8) :: psacws(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: npsacws(mgncol,nlev)    ! number concentration
  ! collection of rain by snow
  real(r8) :: pracs(mgncol,nlev)  ! mass mixing ratio
  real(r8) :: npracs(mgncol,nlev) ! number concentration
  ! freezing of rain
  real(r8) :: mnuccr(mgncol,nlev) ! mass mixing ratio
  real(r8) :: nnuccr(mgncol,nlev) ! number concentration
  ! accretion of droplets by rain
  real(r8) :: pra(mgncol,nlev)    ! mass mixing ratio
  real(r8) :: npra(mgncol,nlev)   ! number concentration
  ! autoconversion of cloud ice to snow
  real(r8) :: prci(mgncol,nlev)   ! mass mixing ratio
  real(r8) :: nprci(mgncol,nlev)  ! number concentration
  ! accretion of cloud ice by snow
  real(r8) :: prai(mgncol,nlev)   ! mass mixing ratio
  real(r8) :: nprai(mgncol,nlev)  ! number concentration
  ! evaporation of rain
  real(r8) :: pre(mgncol,nlev)    ! mass mixing ratio
  ! sublimation of snow
  real(r8) :: prds(mgncol,nlev)   ! mass mixing ratio
  ! number evaporation
  real(r8) :: nsubi(mgncol,nlev)  ! cloud ice
  real(r8) :: nsubc(mgncol,nlev)  ! droplet
  real(r8) :: nsubs(mgncol,nlev)  ! snow
  real(r8) :: nsubr(mgncol,nlev)  ! rain
  ! bergeron process
  real(r8) :: berg(mgncol,nlev)   ! mass mixing ratio (cloud ice)
  real(r8) :: bergs(mgncol,nlev)  ! mass mixing ratio (snow)

  ! fallspeeds
  ! number-weighted
  real(r8) :: uns(mgncol,nlev)    ! snow
  real(r8) :: unr(mgncol,nlev)    ! rain
  ! mass-weighted
  real(r8) :: ums(mgncol,nlev)    ! snow
  real(r8) :: umr(mgncol,nlev)    ! rain
  ! air density corrected fallspeed parameters
  real(r8) :: arn(mgncol,nlev)    ! rain
  real(r8) :: asn(mgncol,nlev)    ! snow
  real(r8) :: acn(mgncol,nlev)    ! cloud droplet
  real(r8) :: ain(mgncol,nlev)    ! cloud ice

  ! saturation vapor pressures
  real(r8) :: esl(mgncol,nlev)    ! liquid
  real(r8) :: esi(mgncol,nlev)    ! ice
  real(r8) :: esn               ! checking for RH after rain evap

  ! saturation vapor mixing ratios
  real(r8) :: qvl(mgncol,nlev)        ! liquid
  real(r8) :: qvi(mgncol,nlev)        ! ice
  real(r8) :: qsn                   ! checking for RH after rain evap

  ! relative humidity
  real(r8) :: relhum(mgncol,nlev)

  ! parameters for cloud water and cloud ice sedimentation calculations
  real(r8) :: fc(nlev)
  real(r8) :: fnc(nlev)
  real(r8) :: fi(nlev)
  real(r8) :: fni(nlev)

  real(r8) :: faloutc(nlev)
  real(r8) :: faloutnc(nlev)
  real(r8) :: falouti(nlev)
  real(r8) :: faloutni(nlev)

  real(r8) :: faltndc
  real(r8) :: faltndnc
  real(r8) :: faltndi
  real(r8) :: faltndni
  real(r8) :: faltndqie
  real(r8) :: faltndqce

  ! sum of source/sink terms for diagnostic precip
  real(r8) :: qstend(mgncol,nlev)    ! snow mixing ratio
  real(r8) :: nstend(mgncol,nlev)     ! snow number concentration
  real(r8) :: qrtend(mgncol,nlev)     ! rain mixing ratio
  real(r8) :: nrtend(mgncol,nlev)     ! rain number concentration
  ! vertically integrated source/sink terms
  real(r8) :: qrtot(mgncol)           ! rain mixing ratio
  real(r8) :: nrtot(mgncol)           ! rain number concentration
  real(r8) :: qstot(mgncol)           ! snow mixing ratio
  real(r8) :: nstot(mgncol)           ! snow number concentration

  ! for calculation of rate1ord_cw2pr_st
  real(r8) :: qcsinksum_rate1ord(mgncol,nlev) ! sum over iterations of cw to precip sink
  real(r8) :: qcsum_rate1ord(mgncol,nlev)     ! sum over iterations of cloud water

  real(r8) :: rainrt(mgncol,nlev)    ! rain rate for reflectivity calculation

  ! dummy variables
  real(r8) :: dum
  real(r8) :: dum1
  ! dummies for checking RH
  real(r8) :: qtmp
  real(r8) :: ttmp
  ! dummies for conservation check
  real(r8) :: qce                   ! qc
  real(r8) :: qie                   ! qi
  real(r8) :: nce                   ! nc
  real(r8) :: nie                   ! ni
  real(r8) :: ratio
  ! dummies for in-cloud variables
  real(r8) :: dumc(mgncol,nlev)   ! qc
  real(r8) :: dumnc(mgncol,nlev)  ! nc
  real(r8) :: dumi(mgncol,nlev)   ! qi
  real(r8) :: dumni(mgncol,nlev)  ! ni
  real(r8) :: dumr(mgncol,nlev)   ! rain mixing ratio
  real(r8) :: dumnr(mgncol,nlev)  ! rain number concentration
  ! Array dummy variable
  real(r8) :: dum_2D(mgncol,nlev)
#ifdef GFDL_COMPATIBLE_MICROP
  real(r8) :: dum2
  real(r8) :: dum2i(mgncol,nlev)   ! used with ice nuleation
  real(r8) :: dum2l(mgncol,nlev)   ! used with drop nuleation
  real(r8) :: dum2a(mgncol,nlev)   ! used with ice nuleation
#endif

  ! loop array variables
  ! "i" and "k" are column/level iterators for internal (MG) variables
  ! "ii" and "kk" are used for indices into input/output buffers
  ! "it" is substepping variable
  ! "n" is used for other iterations (currently just sedimentation)
  integer i, ii, k, kk, it, n

  ! number of iterations for loops over "n"
  integer nstep

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  ! default return error message
  errstring = ' '

  ! Process inputs

  ! assign variable deltat for sub-stepping...
  deltat  = deltatin



#ifdef GFDL_COMPATIBLE_MICROP 
   call pack_diag_array (diag_4d(:,j,:,0:), mgcols, top_lev,   &
                                                       diag_4l(:,j,:,0:))
  call pack_array(    delta_cfin, mgcols, top_lev,   delta_cf )
  call pack_array(   D_eros_l4in, mgcols, top_lev, D_eros_l4  )
  call pack_array(   nerosc4in  , mgcols, top_lev,   nerosc4  )
  call pack_array(   D_eros_i4in, mgcols, top_lev, D_eros_i4  )
  call pack_array(   nerosi4in  , mgcols, top_lev,   nerosi4  )
  call pack_array(    dqcdtin   , mgcols, top_lev,   dqcdt    )
  call pack_array(    dqidtin   , mgcols, top_lev,   dqidt    )
#endif
  call pack_array(          qn, mgcols, top_lev,           q)
  call pack_array(          tn, mgcols, top_lev,           t)
  call pack_array(         qcn, mgcols, top_lev,          qc)
  call pack_array(         qin, mgcols, top_lev,          qi)
  call pack_array(         ncn, mgcols, top_lev,          nc)
  call pack_array(         nin, mgcols, top_lev,          ni)
  call pack_array(          pn, mgcols, top_lev,           p)
  call pack_array(       pdeln, mgcols, top_lev,        pdel)
  call pack_array(     relvarn, mgcols, top_lev,      relvar)
  call pack_array(accre_enhann, mgcols, top_lev, accre_enhan)

  call pack_array(  naain, mgcols, top_lev,  naai)
  call pack_array(npccnin, mgcols, top_lev, npccn)

  ! These are allocated instead of used as automatic arrays
  ! purely to work around a PGI bug.
  allocate(rndst(mgncol,nlev,size(rndstn,3)))
  allocate(nacon(mgncol,nlev,size(rndstn,3)))
  call pack_array( rndstn, mgcols, top_lev, rndst)
  call pack_array(naconin, mgcols, top_lev, nacon)

  if (.not. do_cldice) then
     call pack_array(tnd_qsnown, mgcols, top_lev, tnd_qsnow)
     call pack_array(tnd_nsnown, mgcols, top_lev, tnd_nsnow)
     call pack_array(   re_icen, mgcols, top_lev,    re_ice)
  end if

  ! Some inputs are only used once, to create a reference array that is
  ! used repeatedly later. Rather than bothering to pack these, just
  ! set the local reference directly from the inputs.

  ! pint: used to set int_to_mid
  ! interface to mid-level linear interpolation
  do k = 1,nlev
     do i = 1, mgncol
        ! Set ii and kk to values that correspond to i and k.
        ii = mgcols(i)
        kk = k + top_lev - 1

        int_to_mid(i,k) = (p(i,k) - pint(ii,kk))/ &
             (pint(ii,kk+1) - pint(ii,kk))
     end do
  end do

#ifdef GFDL_COMPATIBLE_MICROP
!initialize gfdl-only arrays-- this may not be needed and will be reviewed later.
    preo(1:mgncol,1:nlev) =0._r8
    prdso(1:mgncol,1:nlev)=0._r8
    cmelo(1:mgncol,1:nlev) =0._r8
    eroslo(1:mgncol,1:nlev) =0._r8
    erosio(1:mgncol,1:nlev) =0._r8
!droplet number
    nucclimo(1:mgncol,1:nlev)   = 0._r8
    npccno(1:mgncol,1:nlev)     = 0._r8
    nnuccco(1:mgncol,1:nlev)    = 0._r8
    nnuccto(1:mgncol,1:nlev)    = 0._r8
    npsacwso(1:mgncol,1:nlev)   = 0._r8
    nsubco(1:mgncol,1:nlev)     = 0._r8
    nerosco(1:mgncol,1:nlev)    = 0._r8
    nprao(1:mgncol,1:nlev)      = 0._r8
    nprc1o(1:mgncol,1:nlev)     = 0._r8
!ice number
    nucclim1io(1:mgncol,1:nlev) = 0._r8
    nnuccdo(1:mgncol,1:nlev)    = 0._r8
    nsacwio(1:mgncol,1:nlev)    = 0._r8
    nsubio(1:mgncol,1:nlev)     = 0._r8
    nerosio(1:mgncol,1:nlev)    = 0._r8
    nprcio(1:mgncol,1:nlev)     = 0._r8
    npraio(1:mgncol,1:nlev)     = 0._r8
#endif

#ifdef GFDL_COMPATIBLE_MICROP
!  initialize bergeron fraction arrays
        sum_freeze(1:mgncol,1:nlev) = 0._r8
        sum_freeze2(1:mgncol,1:nlev) = 0._r8
        sum_rime(1:mgncol,1:nlev)   = 0._r8
        sum_berg(1:mgncol,1:nlev)  = 0._r8
        sum_ice_adj(1:mgncol,1:nlev)  = 0._r8
        sum_bergs(1:mgncol,1:nlev)  = 0._r8
        sum_cond (1:mgncol,1:nlev)  = 0._r8
        sum_splinter(1:mgncol,1:nlev)  = 0._r8
#endif
  ! cldn: used to set cldm, unused for subcolumns
  ! liqcldf: used to set lcldm, unused for subcolumns
  ! icecldf: used to set icldm, unused for subcolumns

  if (microp_uniform) then
     ! subcolumns, set cloud fraction variables to one
     ! if cloud water or ice is present, if not present
     ! set to mincld (mincld used instead of zero, to prevent
     ! possible division by zero errors).

     where (qc >= qsmall)
        lcldm = 1._r8
     elsewhere
        lcldm = mincld
     end where

     where (qi >= qsmall)
        icldm = 1._r8
     elsewhere
        icldm = mincld
     end where

     cldm = max(icldm, lcldm)

  else
     ! get cloud fraction, check for minimum

     do k = 1,nlev
        do i = 1, mgncol
           ii = mgcols(i)
           kk = k + top_lev - 1

           cldm(i,k) = max(cldn(ii,kk),mincld)
           lcldm(i,k) = max(liqcldf(ii,kk),mincld)
           icldm(i,k) = max(icecldf(ii,kk),mincld)
        end do
     end do
  end if

  ! Initialize local variables

  ! local physical properties
  rho = p/(r*t)
  dv = 8.794E-5_r8 * t**1.81_r8 / p
  mu = 1.496E-6_r8 * t**1.5_r8 / (t + 120._r8)
  sc = mu/(rho*dv)

  ! get dz from dp and hydrostatic approx
  ! keep dz positive (define as layer k-1 - layer k)
  dz = pdel/(rho*g)

  ! air density adjustment for fallspeed parameters
  ! includes air density correction factor to the
  ! power of 0.54 following Heymsfield and Bansemer 2007

  rhof=(rhosu/rho)**0.54_r8
#ifdef GFDL_COMPATIBLE_MICROP
  if (.not. rho_factor_in_max_vt) rhof = 1.0
  rhof = MIN (rhof, max_rho_factor_in_vt)
#endif

! --->h1g, add namelist variables, 2014-07-01
! Zhao et al., ACP 2013, Table 1, 
! ai: 350-1400 (s^-1);     as: 5.86-23.44 (m^0.59 s^-1)  
! IceFallFac: 0.5 -- 2;    SnowFallFac: 0.5 -- 2
  arn=ar*rhof
  asn=as*rhof * SnowFallFac

  if( CldFallNew ) then
    acn=g*rhow/(18._r8*mu)
    ain=ai * IceFallFac *(rhosu/rho)**0.35_r8
  else
    acn = ac * rhof
    ain = ai * rhof * IceFallFac
  endif
! <---h1g, 2014-07-01

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! Get humidity and saturation vapor pressures

  do k=1,nlev
     do i=1,mgncol

#ifdef GFDL_COMPATIBLE_MICROP
        call compute_qs (t(i,k), p(i,k), qvl(i,k), q = q(i,k),   &
                             esat = esl(i,k), es_over_liq = .true.)
#else
        call qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k))
#endif

#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (t(i,k), p(i,k), qvi(i,k), q = q(i,k),   &
                             esat = esi(i,k), es_over_liq_and_ice = .true.)
#else
        ! hm fix, make sure when above freezing that esi=esl, not active yet
        if (t(i,k) >= tmelt) then
           esi(i,k)=esl(i,k)
           qvi(i,k)=qvl(i,k)
        else
           call qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k))
        end if
#endif

     end do
  end do

  where (qvl <= 0.0_r8)
     relhum = q
  elsewhere
     relhum = q / min(1.0_r8,qvl)
  end where

  !===============================================
  ! Processes done before substepping
  !===============================================

  ! Initial deposition/sublimation of ice
  !===========================================

  if (do_cldice) then

#ifdef GFDL_COMPATIBLE_MICROP
  do k=1,nlev
   do i=1,mgncol
! define particles available on this step (dum2a).
    if ( liu_in ) then
       dum2a(i,k) = naai(i,k)

! --->h1g, 2014-05-30 add Meyers ice nucleation formula (Only temperature d
! ependent)
    elseif ( use_Meyers ) then
       dum2a(i,k) =  (exp(12.96* (tmelt -t(i,k))/80 - 0.639)) *1000._r8
       dum2a(i,k)= ( dum2a(i,k) )/rho(i,k) ! convert from m-3 to kg-1
! <--- h1g,  2014-05-30    

    elseif (do_ice_nucl_wpdf ) THEN
      if (total_activation) then
        if (activate_all_ice_always) then
          dum2a(i,k) = naai(i,k)
        else
          if (delta_cf(i,k) .gt. 0._r8) then
            dum2a(i,k) = naai(i,k)
          else
            dum2a(i,k)  = 0.
          endif
        endif
      else if (dqa_activation) then
        dum2a(i,k) = (naai(i,k)*max(delta_cf(i,k),0.0) + ni(i,k))/icldm(i,k)
      endif
    else
! default when aerosol field not used for nuclei source
! cooper curve (factor of 1000 is to convert from L-1 to m-3)
      dum2a(i,k)=0.005_r8*exp(0.304_r8*(tmelt -t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
      dum2a(i,k)=min(dum2a(i,k),5.0_r8*exp(0.304_r8*(-tc_cooper)) )/rho(i,k) ! convert from m-3 to kg-1
    endif
  end do
  end do
    call ice_deposition_sublimation_init( clubb_active, deltat, t, q, qc, qi, ni, &
                                          lcldm, icldm, dum2a, rho, dv, &
                                          esl, esi, qvl, qvi, relhum, &
                                       berg, cmei, cmel, &
                                          tiedtke_macrophysics, & 
                                      dqcdt, dqidt)
#else
     
     call ice_deposition_sublimation_init(deltat, t, q, qc, qi, ni, &
                                          lcldm, icldm, naai, rho, dv, &
                                          esl, esi, qvl, qvi, relhum, &
                                          berg, cmei)

#endif
  else
     berg = 0._r8
     cmei = 0._r8
  end if  ! end do_cldice

  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! droplet activation
  ! hm, modify 5/12/11 
  ! get provisional droplet number after activation. This is used for
  ! all microphysical process calculations, for consistency with update of
  ! droplet mass before microphysics

  ! calculate potential for droplet activation if cloud water is present
  ! tendency from activation (npccn) is read in from companion routine
  
  ! output activated liquid and ice (convert from #/kg -> #/m3)
  !--------------------------------------------------
  where (qc >= qsmall)
#ifndef GFDL_COMPATIBLE_MICROP
     nc = nc + npccn*deltat
#endif
     ncal = max(nc/lcldm,cdnl/rho)*rho ! sghan minimum in #/cm3
  elsewhere
     ncal = 0._r8
  end where

  where (t < icenuct)
     ncai = naai*rho
  elsewhere
     ncai = 0._r8
  end where

#ifdef GFDL_COMPATIBLE_MICROP

! calculate ice nucleation dum2i

  do k=1,nlev
   do i=1,mgncol
          if (t(i,k).lt. icenuct       ) then

            if ( liu_in) then
               dum2i(i,k) = naai(i,k)
! --->h1g, 2014-05-30 add Meyers ice nucleation formula (Only temperature dependent)
            elseif ( use_Meyers ) then
               dum2i(i,k) =  (exp(12.96* (tmelt -t(i,k))/80. - 0.639)) *1000._r8
               dum2i(i,k)= ( dum2i(i,k) )/rho(i,k) ! convert from m-3 to kg-1
! <--- h1g,  2014-05-30
            elseif (do_ice_nucl_wpdf ) THEN
 
! if aerosols interact with ice set number of activated ice nuclei
               if (total_activation) then
                 dum2i(i,k) = naai(i,k)
               else if (dqa_activation) then
                 dum2i(i,k) = naai(i,k)
               endif
             else
! cooper curve (factor of 1000 is to convert from L-1 to m-3)
              dum2i(i,k)=0.005_r8*exp(0.304_r8*(tmelt    -t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
            dum2i(i,k)=min(dum2i(i,k),5.0_r8*exp(0.304_r8*(-tc_cooper)))/rho(i,k) ! convert from
!                                                            ! m-3 to kg-1
            endif
          else
            dum2i(i,k)=0._r8
         end if  ! t(i,k).lt. icenuct 
     end do
     end do
#endif

  !INITIALIZE STUFF FOR SUBSTEPPING
  !===============================================

  ! get sub-step time step
  deltat=deltat/real(iter)

  ! hm, set mtime here to avoid answer-changing
  mtime=deltat

  ! initialize tendencies to zero
  tlat1 = 0._r8
  qvlat1 = 0._r8
  qctend1 = 0._r8
  qitend1 = 0._r8
  nctend1 = 0._r8
  nitend1 = 0._r8


  ! initialize microphysics output
  qcsevap=0._r8
  qisevap=0._r8
  qvres  =0._r8
  cmeitot =0._r8
  vtrmc =0._r8
  vtrmi =0._r8
  qcsedten =0._r8
  qisedten =0._r8

  pratot=0._r8
  prctot=0._r8
  mnuccctot=0._r8
  mnuccttot=0._r8
  msacwitot=0._r8
  psacwstot=0._r8
  bergstot=0._r8
  bergtot=0._r8
  melttot=0._r8
  homotot=0._r8
  qcrestot=0._r8
  prcitot=0._r8
  praitot=0._r8
  qirestot=0._r8
  mnuccrtot=0._r8
  pracstot=0._r8
  meltsdttot=0._r8
  frzrdttot=0._r8
  mnuccdtot=0._r8

  rflx=0._r8
  sflx=0._r8

  ! initialize precip output
  
  qrout=0._r8
  qsout=0._r8
  nrout=0._r8
  nsout=0._r8

  ! for refl calc
  rainrt = 0._r8
  
#ifdef GFDL_COMPATIBLE_MICROP
  dum2l = 0.
#endif

  ! initialize rain size
  rercld=0._r8
  arcld = 0

  qcsinksum_rate1ord = 0._r8 
  qcsum_rate1ord     = 0._r8 

  ! initialize variables for trop_mozart
  nevapr = 0._r8
  evapsnow = 0._r8
  prain = 0._r8
  prodsnow = 0._r8
  cmeout = 0._r8

  prect1 = 0._r8
  preci1 = 0._r8

  cldmax = mincld

  lamc=0._r8

#ifdef GFDL_COMPATIBLE_MICROP
  cmel_orig = cmel
  cmei_orig = cmei
  berg_orig = berg
#endif

  !*********DO SUBSTEPPING!***************
  !============================================
  substepping: do it=1,iter

     ! initialize sub-step microphysical tendencies

     tlat=0._r8
     qvlat=0._r8
     qctend=0._r8
     qitend=0._r8
     qstend = 0._r8
     qrtend = 0._r8
     nctend=0._r8
     nitend=0._r8
     nrtend = 0._r8
     nstend = 0._r8

     ! initialize diagnostic precipitation to zero
     qcic  = 0._r8
     qiic  = 0._r8
     qsic  = 0._r8
     qric  = 0._r8

     ncic  = 0._r8
     niic  = 0._r8
     nsic  = 0._r8
     nric  = 0._r8

     ! initialize precip at surface

     prect = 0._r8
     preci = 0._r8

     ! initialize vertically-integrated rain and snow tendencies

     qrtot = 0._r8
     nrtot = 0._r8
     qstot = 0._r8
     nstot = 0._r8

     ! recalculate saturation vapor pressure for liquid and ice
     do k = 1, nlev
        do i = 1, mgncol

#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (t(i,k), p(i,k), qvl(i,k), q = q(i,k),   &
                             esat = esl(i,k), es_over_liq = .true.)
#else
           call qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k))
#endif

#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (t(i,k), p(i,k), qvi(i,k), q = q(i,k),   &
                             esat = esi(i,k), es_over_liq_and_ice = .true.)
#else
           ! hm fix, make sure when above freezing that esi=esl, not active yet
           if (t(i,k) >= tmelt) then
              esi(i,k)=esl(i,k)
              qvi(i,k)=qvl(i,k)
           else
              call qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k))
           end if
#endif

        end do
     end do

     where (qvl <= 0.0_r8)
        relhum = q
     elsewhere
        relhum = q / min(1.0_r8,qvl)
     end where

#ifdef GFDL_COMPATIBLE_MICROP
!  set / return  erosion, bergeron and condensation fields to input values
          if( .not. clubb_active) then
            nerosi = nerosi4
            nerosc = nerosc4
            D_eros_l = D_eros_l4
            D_eros_i = D_eros_i4
            cmel = cmel_orig
            cmei = cmei_orig
            berg = berg_orig
          endif
#endif

     ! decrease in number concentration due to sublimation/evap
     !-------------------------------------------------------
     ! divide by cloud fraction to get in-cloud decrease
     ! don't reduce Nc due to bergeron process

#ifdef GFDL_COMPATIBLE_MICROP
      if ( .not. tiedtke_macrophysics) then
        where (cmei < 0._r8 .and. qi > qsmall .and. icldm > mincld)
          nsubi = cmei / qi * ni / icldm
        elsewhere
          nsubi = 0._r8
        end where
      else
          nsubi = 0._r8
      endif

#else
     where (cmei < 0._r8 .and. qi > qsmall .and. icldm > mincld)
        nsubi = cmei / qi * ni / icldm
     elsewhere
        nsubi = 0._r8
     end where
#endif

#ifdef GFDL_COMPATIBLE_MICROP
      if ( .not. tiedtke_macrophysics) then
!!SHOULD NSUBC be zero ?? NCAR says yes, MG includes code but it is 
!! not activated with Tiedtke. Should it be activated in general for 
!! non-Tiedtke case ??
        if( clubb_active ) then
           nsubc = 0._r8
        else  
        where (cmel < 0._r8 .and. qc > qsmall .and. cldm > mincld)
           nsubc = cmel / qc * nc / cldm
        elsewhere
           nsubc = 0._r8
        end where
        endif   !(clubb_active)
     else  ! tiedtke_macrophysics
           nsubc = 0._r8
     endif
#else
     nsubc = 0._r8
#endif

#ifdef GFDL_COMPATIBLE_MICROP
!  Note that ice nucleation calculated later on in code for 
!  Tiedtke macrophysics case.


if ( ActNew ) then

  if (.not. tiedtke_macrophysics) then
     if (do_cldice) then
      
! --->h1g, 2014-05-30 add Meyers ice nucleation option
      if( liu_in .or. use_Meyers ) then
! <--- h1g,  2014-05-30
        where (dum2i > 0._r8 .and. t < icenuct .and. &
          relhum*esl/esi > rhmini+0.05_r8)
                  nnuccd = (dum2i - ni/icldm)/deltat*icldm
                  nnuccd = max(nnuccd, 0._r8)
        elsewhere
           nnuccd = 0._r8
           nimax = 0._r8
           mnuccd = 0._r8
        end where
      elseif (total_activation) then
        where (dum2i > 0._r8 .and. t < icenuct .and. &
          relhum*esl/esi > rhmini+0.05_r8)

           !if NAAI > 0. then set numice = naai (as before)
           !note: this is gridbox averaged
           ! hm, modify to use mtime
           nnuccd = (dum2i-ni/icldm)/mtime*icldm
           nnuccd = max(nnuccd,0._r8)
        elsewhere
           nnuccd = 0._r8
           nimax = 0._r8
           mnuccd = 0._r8
        end where


     else if (dqa_activation) then
        where (dum2i > 0._r8 .and. t < icenuct .and. &
          relhum*esl/esi > rhmini+0.05_r8)

           !if NAAI > 0. then set numice = naai (as before)
           !note: this is gridbox averaged
           ! hm, modify to use mtime
           nnuccd = max(delta_cf,0._r8) *dum2i/deltatin
        elsewhere
           nnuccd = 0._r8
           nimax = 0._r8
           mnuccd = 0._r8
        end where
!---> h1g, for cooper nucleation, 2014-05-29
      else
        where (dum2i > 0._r8 .and. t < icenuct .and. &
          relhum*esl/esi > rhmini+0.05_r8)
                  nnuccd = (dum2i - ni/icldm)/deltat*icldm
                  nnuccd = max(nnuccd, 0._r8)
        elsewhere
           nnuccd = 0._r8
           nimax = 0._r8
           mnuccd = 0._r8
        end where
!<--- h1g, 2014-05-29
      
      endif ! liu_in .or. use_Meyers
           nimax = dum2i*icldm

         if (.not. dqa_activation) then
           !Calc mass of new particles using new crystal mass...
           !also this will be multiplied by mtime as nnuccd is...

           mnuccd = nnuccd * mi0

           !  add mnuccd to cmei....
           cmei = cmei + mnuccd
        
           !  limit cmei
           !-------------------------------------------------------
           cmei = min(cmei,(q-qvi)/calc_ab(t, qvi, xxls)/deltat)

           ! limit for roundoff error
           cmei = cmei * omsm
         else
           mnuccd = 0._r8
         endif

     end if  ! do_cldice
     end if  ! .not. tiedtke
 else !use cldwat2m_micro.F90 activation
  do k = 1, nlev
    do i = 1, mgncol

      if( .not. tiedtke_macrophysics) then
        if (do_cldice) then
          if (dum2i(i,k).gt.0._r8.and.t(i,k).lt. icenuct       .and. &
                    relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then
            if( liu_in .or. use_Meyers ) then
              nnuccd(i,k) = (dum2i(i,k) - ni(i,k)/icldm(i,k))/deltat*  &
                                                                icldm(i,k)
              nnuccd(i,k) = max(nnuccd(i,k), 0._r8)
            elseif (total_activation) then
              nnuccd(i,k) = (dum2i(i,k) - ni(i,k)/icldm(i,k))/mtime*  &
                                                                icldm(i,k)
              nnuccd(i,k) = max(nnuccd(i,k), 0._r8)
            else if (dqa_activation) then
              nnuccd(i,k) = max(delta_cf(i,k),0._r8) *dum2i(i,k)/deltatin
            else  ! for cooper nucleation, 2014-05-29
              nnuccd(i,k) = (dum2i(i,k) - ni(i,k)/icldm(i,k))/deltat*  &
                                                                icldm(i,k)
              nnuccd(i,k) = max(nnuccd(i,k), 0._r8)
            endif !  liu_in .or. use_Meyers
            nimax(i,k) = dum2i(i,k)*icldm(i,k)

            if (.not. dqa_activation) then
            !Calc mass of new particles using new crystal mass...
            !also this will be multiplied by mtime as nnuccd is...
              mnuccd(i,k) = nnuccd(i,k) * mi0
            !  add mnuccd to cmei....
              cmei(i,k) = cmei(i,k) + mnuccd(i,k)
            !  limit cmei
            !-------------------------------------------------------
              cmei(i,k) = min(cmei(i,k),(q(i,k)-qvi(i,k))/calc_ab(t(i,k), qvi(i,k), xxls)/deltat)

            ! limit for roundoff error
              cmei(i,k) = cmei(i,k) * omsm
            else
              mnuccd(i,k) = 0._r8
            endif ! dqa_activation
 
          else ! NOT (  dum2i(i,k).gt.0._r8.and.t(i,k).lt. icenuct  ... )                         
            nnuccd(i,k) = 0._r8
            nimax(i,k)  = 0._r8
            mnuccd(i,k) = 0._r8
         end if  ! (  dum2i(i,k).gt.0._r8.and.t(i,k).lt. icenuct  ... )         
        end if  ! do_cldice
      end if  ! .not. tiedtke

    enddo ! do i = 1, mgncol
  enddo ! do k = 1, nlev

endif !ActNew

#else
     ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
     !-------------------------------------------------------

     if (do_cldice) then
        where (naai > 0._r8 .and. t < icenuct .and. &
          relhum*esl/esi > rhmini+0.05_r8)

           !if NAAI > 0. then set numice = naai (as before)
           !note: this is gridbox averaged
           ! hm, modify to use mtime
           nnuccd = (naai-ni/icldm)/mtime*icldm
           nnuccd = max(nnuccd,0._r8)
           nimax = naai*icldm

           !Calc mass of new particles using new crystal mass...
           !also this will be multiplied by mtime as nnuccd is...

           mnuccd = nnuccd * mi0

           !  add mnuccd to cmei....
           cmei = cmei + mnuccd
        
           !  limit cmei
           !-------------------------------------------------------
           cmei = min(cmei,(q-qvi)/calc_ab(t, qvi, xxls)/deltat)

           ! limit for roundoff error
           cmei = cmei * omsm

        elsewhere
           nnuccd = 0._r8
           nimax = 0._r8
           mnuccd = 0._r8
        end where

     end if

#endif
     pre_vert_loop: do k=1,nlev

        pre_col_loop: do i=1,mgncol

           ! for microphysical process calculations
           ! units are kg/kg for mixing ratio, 1/kg for number conc

#ifdef GFDL_COMPATIBLE_MICROP
!  adjust previously calculated tendencies to avoid creating negative
!  water species
     
         if ( .not. clubb_active) then
            if (qc(i,k) + (cmel(i,k) + D_eros_l(i,k) -    &
                                         berg(i,k))*deltat.lt.qsmall) then
              qcic(i,k)=0._r8
              ncic(i,k)=0._r8
              if (qc(i,k) + (cmel(i,k) + D_eros_l(i,k) -   &
                                         berg(i,k))*deltat.lt.0._r8) then
                if (cmel(i,k).lt.0._r8) then
!++ first only scale cmel, d_eros
                  dum = -cmel(i,k) - D_eros_l(i,k)
                  if (dum .gt. 1.e-30_r8) then
                    dum3 = qc(i,k)/deltat/dum*omsm
                  else
                    dum3 = 0._r8
                  end if
                  cmel(i,k) = dum3*cmel(i,k)
                  D_eros_l(i,k) = dum3*D_eros_l(i,k)
                  dum = -cmel(i,k) - D_eros_l(i,k) + berg(i,k)
                  if (dum .gt. 1.e-30_r8) then
                    dum3 = qc(i,k)/deltat/dum*omsm
                  else
                    dum3 = 0._r8
                  end if
                  cmel(i,k) = dum3*cmel(i,k)
                  D_eros_l(i,k) = dum3*D_eros_l(i,k)
                  berg(i,k) = dum3*berg(i,k)
                else
                  dum = -D_eros_l(i,k) + berg(i,k)
                  if (dum .gt. 1.e-30_r8) then
                    dum3 = ( qc(i,k)/deltat +  cmel(i,k) ) / dum * omsm
                  else
                    dum3 = 0._r8
                  end if
                  D_eros_l(i,k) = D_eros_l(i,k)*dum3
                  berg(i,k) = berg(i,k)*dum3
                endif
              endif
            else
              qcic(i,k)=min(qc(i,k)/lcldm(i,k),5.e-3_r8)
              ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
!-->cjg
! hm add 6/2/11 specify droplet concentration
               if (nccons) then
               ncic(i,k)=ncnst/rho(i,k)
               end if
            end if
           
            if (qi(i,k) + (cmei(i,k) + D_eros_i(i,k) + berg(i,k))*   &
                                                     deltat.lt.qsmall) then
              qiic(i,k)=0._r8
              niic(i,k)=0._r8
              if (qi(i,k) + (cmei(i,k) + D_eros_i(i,k) + berg(i,k))*  &
                                                     deltat.lt.0._r8) then
                if (cmei(i,k).lt.0._r8) then
                  dum = - cmei(i,k) - D_eros_i(i,k)
                  if (dum .gt. 1.e-30_r8) then
                    dum3 = (qi(i,k)/deltat + berg(i,k))/dum*omsm
                  else
                    dum3 = 0._r8
                  end if
                  cmei(i,k) = dum3 * cmei(i,k)
                  D_eros_i(i,k) = dum3 *  D_eros_i(i,k)
               else
                  dum = - D_eros_i(i,k)
                  if (dum .gt. 1.e-30_r8) then
                    dum3 = (qi(i,k)/deltat + cmei(i,k) + berg(i,k))/  &
                                                                dum*omsm
                  else
                    dum3 = 0._r8
                   end if
                  D_eros_i(i,k) = dum3*D_eros_i(i,k)
                end if
              end if
            else
              qiic(i,k)=min(qi(i,k)/icldm(i,k),5.e-3_r8)
              niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)
!-->cjg
! hm add 6/2/11 switch for specification of cloud ice number
           if (nicons) then
           niic(i,k)=ninst/rho(i,k)
           end if             
            end if
         else
#endif 
           if (qc(i,k) - berg(i,k)*deltat.ge.qsmall) then
              ! limit in-cloud values to 0.005 kg/kg
              qcic(i,k)=min(qc(i,k)/lcldm(i,k),5.e-3_r8)
              ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)

              ! hm add 6/2/11 specify droplet concentration
              if (nccons) then
                 ncic(i,k)=ncnst/rho(i,k)
              end if
           else
              qcic(i,k)=0._r8
              ncic(i,k)=0._r8

              berg(i,k)=qc(i,k)/deltat*omsm
           end if

           if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.ge.qsmall) then
              ! limit in-cloud values to 0.005 kg/kg
              qiic(i,k)=min(qi(i,k)/icldm(i,k),5.e-3_r8)
              niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)

              ! hm add 6/2/11 switch for specification of cloud ice number
              if (nicons) then
                 niic(i,k)=ninst/rho(i,k)
              end if
           else
              qiic(i,k)=0._r8
              niic(i,k)=0._r8

              if (do_cldice) then
                 cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm
              end if
           end if
#ifdef GFDL_COMPATIBLE_MICROP
         endif  ! (clubb_active)
#endif

        end do pre_col_loop
     end do pre_vert_loop

#ifdef GFDL_COMPATIBLE_MICROP
! calculate ice nucleation here for tiedtke case, since it needs qcic 
! which was not defined until the previous  loop.
            if (tiedtke_macrophysics) then

     if (do_cldice) then
      do k=1,nlev

      do i=1,mgncol
           if (qiic(i,k) > qsmall .and. t(i,k) < icenuct ) then 

           !if NAAI > 0. then set numice = naai (as before)
           !note: this is gridbox averaged
           ! hm, modify to use mtime
                if (total_activation) then
                  if (activate_all_ice_always) then
                    nnuccd(i,k) = (dum2i(i,k)-ni(i,k)/icldm(i,k))/mtime*icldm(i,k)
                  else 
                       if (delta_cf(i,k) .gt. 0._r8) then
                       nnuccd(i,k) = (dum2i(i,k) - ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
                     else
                       nnuccd(i,k) = 0._r8
                     endif    
                   endif
                   nnuccd(i,k) = max(nnuccd(i,k),0._r8)
                else if (dqa_activation) then
                  nnuccd(i,k) = max(delta_cf(i,k),0._r8) *dum2i(i,k)/deltatin
                 endif
                nimax(i,k) = dum2i(i,k)*icldm(i,k)


        else
           nnuccd(i,k) = 0._r8
           nimax(i,k) = 0._r8
           mnuccd(i,k) = 0._r8
        end if
     end do
     end do

     end if
     end if
#endif
     ! add to cme output
     cmeout = cmeout + cmei

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! droplet activation

#ifdef GFDL_COMPATIBLE_MICROP
   do k=1,nlev
    do i=1,mgncol
     if (qcic(i,k).ge.qsmall) then
        dum2l(i,k) = max(0._r8, npccn  (i,k))
      if (clubb_active) then
!RSH npccn2 is the change in droplet number on this step. In the 
! non-GFDL_COMPATIBLE_MICROP, this is not calculated since the input 
! droplet number already includes the newly activated droplets.
        npccn2(i,k) = ((dum2l  (i,k) - nc(i,k)/cldm(i,k))/deltat)* &
                                                                  cldm(i,k)
        npccn2(i,k) = max(0._r8,npccn2(i,k))
      else  ! (clubb_active)
        IF ( total_activation) THEN
          npccn2(i,k) = ((dum2l  (i,k) - nc(i,k)/cldm(i,k))/deltat)* &
                                                                  cldm(i,k)
          npccn2(i,k) = max(0._r8,npccn2(i,k))
                                                             ! in #/cm3  
       ELSE IF   ( dqa_activation    ) THEN
!delta_cf:  A_dt * (1.-qabar)   where A_dt = A*dt , A source rate

!dum2l has already been multiplied by 1.e6/airdens(i,k)
         npccn2(i,k) = max (delta_cf(i,k), 0._r8)*dum2l  (i,k)/deltatin
       END IF
        ncmax(i,k) = dum2l  (i,k)*cldm(i,k)
       endif ! (clubb_active)
    else
       npccn2(i,k)=0._r8
       ncmax(i,k) = 0._r8
    end if  
      end do
      end do
#endif

     !=========================================================
     ! Main microphysical loop
     !=========================================================

     ! initialize precip fallspeeds to zero
     ums = 0._r8
     uns = 0._r8
     umr = 0._r8
     unr = 0._r8

     ! for sub-columns cldm has already been set to 1 if cloud
     ! water or ice is present, so cldmax will be correctly set below
     ! and nothing extra needs to be done here

     cldmax = cldm

     micro_vert_loop: do k=1,nlev

        ! calculate precip fraction based on maximum overlap assumption

        ! if rain or snow mix ratios are smaller than threshold, 
        ! then leave cldmax as cloud fraction at current level
        if (k /= 1) then
           where (qric(:,k-1).ge.qsmall .or. qsic(:,k-1).ge.qsmall)
              cldmax(:,k)=max(cldmax(:,k-1),cldmax(:,k))
           end where
        end if
        
        do i = 1, mgncol

           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           ! get size distribution parameters based on in-cloud cloud water
           ! these calculations also ensure consistency between number and mixing ratio
           !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

           ! cloud liquid
           !-------------------------------------------

           ! ".true." below turns on adjustment of ncic for consistency.
           call size_dist_param_liq(qcic(i,k), ncic(i,k), cdnl, rho(i,k), .true., &
                pgam(i,k), lamc(i,k))

           if (lamc(i,k) > 0._r8) then

              ! parameter to calculate droplet freezing
              cdist1(i,k) = ncic(i,k)/gamma(pgam(i,k)+1._r8)

           else
              cdist1(i,k) = 0._r8
           end if

        end do

        !========================================================================
        ! autoconversion of cloud liquid water to rain
        ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
        ! minimum qc of 1 x 10^-8 prevents floating point error

#ifdef GFDL_COMPATIBLE_MICROP
       do i=1,mgncol
        call kk2000_liq_autoconversion(qcic(i,k), ncic(i,k), rho(i,k), &
             relvar(i,k), prc(i,k), nprc(i,k), nprc1(i,k))
       end do
#else
        call kk2000_liq_autoconversion(qcic(:,k), ncic(:,k), rho(:,k), &
             relvar(:,k), prc(:,k), nprc(:,k), nprc1(:,k))
#endif

        ! add autoconversion to precip from above to get provisional rain mixing ratio
        ! and number concentration (qric and nric)

        ! hm 11-16-11, modify, divide dz by 2 to use modified mid-point method
        ! This estimates rain and snow mass and number mixing ratios at
        ! mid-point to calculate process rates at mid-point, with final
        ! values of rain and snow mass and number mixing ratios calculated
        ! on interfaces

        if (k .eq. 1) then
           dum=0.45_r8

           qric(:,k)= prc(:,k)*lcldm(:,k)*dz(:,k)/2._r8/cldmax(:,k)/dum
           nric(:,k)=nprc(:,k)*lcldm(:,k)*dz(:,k)/2._r8/cldmax(:,k)/dum
        else

#ifdef GFDL_COMPATIBLE_MICROP
              if (allow_all_cldtop_collection) then
! NCAR allows no autoconversion of rain number if rain/snow falling from 
! above. this assumes that new drizzle drops formed by autoconversion are 
! rapidly collected by the existing rain/snow particles falling from above.
! Marc's code allowed autoconversion to change rain number, so  variable 
! allow_all_cldtop_collection  introduced, which when .true. would turn off
! this effect. By default, it is .false. for GFDL (as in MG), in contrast 
! to what NCAR does (ifndef GFDL_COMPATIBLE_MICROP). 
             where (qric(:,k-1).ge.1.e-9_r8 .or.    &
                                        qsic(:,k-1).ge.1.e-9_r8) 
                  nprc(:,k)=0._r8
                end where
              endif  !  allow_all_cldtop_collection
#else


           ! no autoconversion of rain number if rain/snow falling from above
           ! this assumes that new drizzle drops formed by autoconversion are rapidly collected
           ! by the existing rain/snow particles from above

           where (qric(:,k-1).ge.1.e-9_r8 .or. qsic(:,k-1).ge.1.e-9_r8)
              nprc(:,k) = 0._r8
           end where
#endif

           do i = 1,mgncol
              if (qric(i,k-1).ge.qsmall) then
                 dum=umr(i,k-1)
                 dum1=unr(i,k-1)
              else
                 ! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)

                 dum=0.45_r8
                 dum1=0.45_r8
              end if

              qric(i,k) = (rho(i,k-1)*umr(i,k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
                   (rho(i,k)*dz(i,k)/2._r8*((pra(i,k-1)+prc(i,k))*lcldm(i,k)+ &
                   (pre(i,k-1)-pracs(i,k-1)-mnuccr(i,k-1))*cldmax(i,k)))) &
                   /(dum*rho(i,k)*cldmax(i,k))
              nric(i,k) = (rho(i,k-1)*unr(i,k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
                   (rho(i,k)*dz(i,k)/2._r8*(nprc(i,k)*lcldm(i,k)+ &
                   (nsubr(i,k-1)-npracs(i,k-1)-nnuccr(i,k-1)+nragg(i,k-1))*cldmax(i,k))))&
                   /(dum1*rho(i,k)*cldmax(i,k))

           end do

        end if

        ! if precip mix ratio is zero so should number concentration

        where (qric(:,k).lt.qsmall)
           qric(:,k)=0._r8
           nric(:,k)=0._r8
        end where

        ! make sure number concentration is a positive number to avoid 
        ! taking root of negative later

        nric(:,k)=max(nric(:,k),0._r8)

        ! Get size distribution parameters for cloud ice

        call size_dist_param_ice(qiic(:,k), niic(:,k), lami(:,k), n0i(:,k))

        !.......................................................................
        ! Autoconversion of cloud ice to snow
        ! similar to Ferrier (1994)

        if (do_cldice) then
           call ice_autoconversion(t(:,k), qiic(:,k), lami(:,k), n0i(:,k), &
                prci(:,k), nprci(:,k))
           
        else
           ! Add in the particles that we have already converted to snow, and
           ! don't do any further autoconversion of ice.
           prci(:,k)  = tnd_qsnow(:,k) / cldm(:,k)
           nprci(:,k) = tnd_nsnow(:,k) / cldm(:,k)
        end if

        do i=1,mgncol

           ! add autoconversion to flux from level above to get provisional snow mixing ratio
           ! and number concentration (qsic and nsic)

           ! hm 11-16-11 modify for mid-point method, see comments above

           if (k == 1) then
              dum=(asn(i,k)*cons25)
              qsic(i,k)=prci(i,k)*icldm(i,k)*dz(i,k)/2._r8/cldmax(i,k)/dum
              nsic(i,k)=nprci(i,k)*icldm(i,k)*dz(i,k)/2._r8/cldmax(i,k)/dum
           else
              if (qsic(i,k-1) >= qsmall) then
                 dum=ums(i,k-1)
                 dum1=uns(i,k-1)
              else
                 dum = asn(i,k)*cons25
                 dum1 = dum
              end if

              qsic(i,k) = (rho(i,k-1)*ums(i,k-1)*qsic(i,k-1)*cldmax(i,k-1)+ &
                   (rho(i,k)*dz(i,k)/2._r8*((prci(i,k)+prai(i,k-1)+psacws(i,k-1)+bergs(i,k-1))*icldm(i,k)+ &
                   (prds(i,k-1)+pracs(i,k-1)+mnuccr(i,k-1))*cldmax(i,k))))&
                   /(dum*rho(i,k)*cldmax(i,k))

              nsic(i,k) = (rho(i,k-1)*uns(i,k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
                   (rho(i,k)*dz(i,k)/2._r8*(nprci(i,k)*icldm(i,k)+ &
                   (nsubs(i,k-1)+nsagg(i,k-1)+nnuccr(i,k-1))*cldmax(i,k)))) &
                   /(dum1*rho(i,k)*cldmax(i,k))

           end if

        end do

        ! if precip mix ratio is zero so should number concentration

        where (qsic(:,k) < qsmall)
           qsic(:,k)=0._r8
           nsic(:,k)=0._r8
        end where

        ! make sure number concentration is a positive number to avoid 
        ! taking root of negative later

        nsic(:,k)=max(nsic(:,k),0._r8)

        !.......................................................................
        ! get size distribution parameters for precip
        !......................................................................
        ! rain

        call size_dist_param_rain(qric(:,k), nric(:,k), lamr(:,k), n0r(:,k))

        where (lamr(:,k) >= qsmall)

           ! provisional rain number and mass weighted mean fallspeed (m/s)

           unr(:,k) = min(arn(:,k)*cons4/lamr(:,k)**br,9.1_r8*rhof(:,k))
           umr(:,k) = min(arn(:,k)*cons5/(6._r8*lamr(:,k)**br),9.1_r8*rhof(:,k))

        elsewhere
           umr(:,k) = 0._r8
           unr(:,k) = 0._r8
        end where

        !......................................................................
        ! snow

        call size_dist_param_snow(qsic(:,k), nsic(:,k), lams(:,k), n0s(:,k))

        where (lams(:,k) > 0._r8)

           ! provisional snow number and mass weighted mean fallspeed (m/s)

           ums(:,k) = min(asn(:,k)*cons8/(6._r8*lams(:,k)**bs),1.2_r8*rhof(:,k))
           uns(:,k) = min(asn(:,k)*cons7/lams(:,k)**bs,1.2_r8*rhof(:,k))

        elsewhere
           ums(:,k) = 0._r8
           uns(:,k) = 0._r8
        end where

        if (do_cldice) then

           ! heterogeneous freezing of cloud water
           !----------------------------------------------

#ifdef GFDL_COMPATIBLE_MICROP
      do i=1,mgncol
           call immersion_freezing(t(i,k), pgam(i,k), lamc(i,k), cdist1(i,k), qcic(i,k), &
                relvar(i,k), mnuccc(i,k), nnuccc(i,k))
       end do
#else
           call immersion_freezing(t(:,k), pgam(:,k), lamc(:,k), cdist1(:,k), qcic(:,k), &
                relvar(:,k), mnuccc(:,k), nnuccc(:,k))
#endif

           ! make sure number of droplets frozen does not exceed available ice nuclei concentration
           ! this prevents 'runaway' droplet freezing


#ifdef GFDL_COMPATIBLE_MICROP
           where (qcic(:,k).ge.qsmall .and. t(:,k).lt.tmelt-4._r8)
#else
           where (qcic(:,k).ge.qsmall .and. t(:,k).lt.269.15_r8)
#endif
              where (nnuccc(:,k)*lcldm(:,k).gt.nnuccd(:,k))
                 ! scale mixing ratio of droplet freezing with limit
                 mnuccc(:,k)=mnuccc(:,k)*(nnuccd(:,k)/(nnuccc(:,k)*lcldm(:,k)))
                 nnuccc(:,k)=nnuccd(:,k)/lcldm(:,k)
              end where
           end where

#ifdef GFDL_COMPATIBLE_MICROP
!   contact freezing not currently availablein GFDL. Need to get proper
!   dust input fields.
             if (.not. clubb_active) then
              mnucct(:,k) = 0._r8
              nnucct(:,k) = 0._r8
             else
           call contact_freezing(t(:,k), p(:,k), rndst(:,k,:), nacon(:,k,:), &
                pgam(:,k), lamc(:,k), cdist1(:,k), qcic(:,k), &
                relvar(:,k), mnucct(:,k), nnucct(:,k))
             endif
#else  
           call contact_freezing(t(:,k), p(:,k), rndst(:,k,:), nacon(:,k,:), &
                pgam(:,k), lamc(:,k), cdist1(:,k), qcic(:,k), &
                relvar(:,k), mnucct(:,k), nnucct(:,k))
#endif

        else
           mnuccc(:,k)=0._r8
           nnuccc(:,k)=0._r8
           mnucct(:,k)=0._r8
           nnucct(:,k)=0._r8
        end if

        call snow_self_aggregation(t(:,k), rho(:,k), asn(:,k), qsic(:,k), nsic(:,k), &
             nsagg(:,k))

        call accrete_cloud_water_snow(t(:,k), rho(:,k), asn(:,k), uns(:,k), mu(:,k), &
             qcic(:,k), ncic(:,k), qsic(:,k), pgam(:,k), lamc(:,k), lams(:,k), n0s(:,k), &
             psacws(:,k), npsacws(:,k))

#ifdef GFDL_COMPATIBLE_MICROP
            if (do_hallet_mossop .or. clubb_active ) then
#endif
        if (do_cldice) then
           call secondary_ice_production(t(:,k), psacws(:,k), msacwi(:,k), nsacwi(:,k))
        else
           nsacwi(:,k) = 0.0_r8
           msacwi(:,k) = 0.0_r8
        end if
#ifdef GFDL_COMPATIBLE_MICROP
            else
              msacwi(:,k) = 0.0_r8
              nsacwi(:,k) = 0.0_r8
            endif ! (do_hallet_mossop .or. clubb_active)
#endif

        call accrete_rain_snow(t(:,k), rho(:,k), umr(:,k), ums(:,k), unr(:,k), uns(:,k), &
             qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
             pracs(:,k), npracs(:,k))

        call heterogeneous_rain_freezing(t(:,k), qric(:,k), nric(:,k), lamr(:,k), &
             mnuccr(:,k), nnuccr(:,k))

#ifdef GFDL_COMPATIBLE_MICROP
       if(clubb_active) then
        do i=1,mgncol
           if( use_qcvar_in_accretion ) then
             if ( relvar(i,k) > qcvar_max4accr ) then
                 accr_scale = 1.0
             elseif( relvar(i,k) < qcvar_min4accr ) then
                 accr_scale = accretion_scale_max
             else
                 accr_scale =   (accretion_scale_max-1.0)/  &
                  (1.0/qcvar_min4accr - 1.0/qcvar_max4accr) &
                * (1.0/relvar(i,k) - 1.0/qcvar_max4accr     ) + 1.0
             endif
           else
             accr_scale = accretion_scale
           endif
        call accrete_cloud_water_rain(qric(i,k), qcic(i,k), ncic(i,k), &
            rho(i,k), umr(i,k), relvar(i,k), accr_scale, pra(i,k),   &
            npra(i,k))
        end do
       else
        do i=1,mgncol
! ---> h1g, boost accretion in MG, 2014-07-17
         if( use_const_accretion_scale ) then
           call accrete_cloud_water_rain(qric(i,k), qcic(i,k), ncic(i,k), &
              rho(i,k), umr(i,k), relvar(i,k), accretion_scale, pra(i,k), &
                                                                npra(i,k))
         else
           call accrete_cloud_water_rain(qric(i,k), qcic(i,k), ncic(i,k), &
              rho(i,k), umr(i,k), relvar(i,k), accre_enhan(i,k),    &
                                                      pra(i,k), npra(i,k))
         endif
! <--- h1g, 2014-07-17
        end do
       endif
#else
        call accrete_cloud_water_rain(qric(:,k), qcic(:,k), ncic(:,k), &
             rho(i,k), umr(i,k), &
             relvar(:,k), accre_enhan(:,k), pra(:,k), npra(:,k))
#endif

        call self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k))

        if (do_cldice) then
           call accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), &
                qsic(:,k), lams(:,k), n0s(:,k), prai(:,k), nprai(:,k))
        else
           prai(:,k) = 0._r8
           nprai(:,k) = 0._r8
        end if
              
#ifdef GFDL_COMPATIBLE_MICROP
        do i=1,mgncol
        call evaporate_sublimate_precip(deltat, t(i,k), p(i,k), rho(i,k), &
             dv(i,k), mu(i,k), sc(i,k), q(i,k), qvl(i,k), qvi(i,k), &
             lcldm(i,k), cldmax(i,k), arn(i,k), asn(i,k), qcic(i,k), qiic(i,k), &
             qric(i,k), qsic(i,k), lamr(i,k), n0r(i,k), lams(i,k), n0s(i,k), cmei(i,k), &
             D_eros_l(i,k), D_eros_i(i,k), cmel(i,k), &
             pre(i,k), prds(i,k),clubb_active)
         end do
#else
        call evaporate_sublimate_precip(deltat, t(:,k), p(:,k), rho(:,k), &
             dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
             lcldm(:,k), cldmax(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
             qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), cmei(:,k), &
             pre(:,k), prds(:,k))
#endif

        call bergeron_process(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), &
             qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), qsic(:,k), lams(:,k), n0s(:,k), &
             bergs(:,k))

        ! Big "administration" loop enforces conservation, updates variables
        ! that accumulate over substeps, and sets output variables.

        do i=1,mgncol

           ! conservation to ensure no negative values of cloud water/precipitation
           ! in case microphysical process rates are large
           !===================================================================

           ! make sure to use end-of-time step values for cloud water, ice, due
           ! condensation/deposition

           ! note: for check on conservation, processes are multiplied by omsm
           ! to prevent problems due to round off error

#ifdef GFDL_COMPATIBLE_MICROP
        if (.not. clubb_active) then
           qce = (qc(i,k) + (D_eros_l(i,k) + cmel(i,k) - berg(i,k))*deltat)
           qie = (qi(i,k) + (D_eros_i(i,k) + cmei(i,k) + berg(i,k))*deltat)
        else
#endif
           qce = (qc(i,k) - berg(i,k)*deltat)
           qie = (qi(i,k) + (cmei(i,k) + berg(i,k))*deltat)
#ifdef GFDL_COMPATIBLE_MICROP
        endif
#endif
#ifdef GFDL_COMPATIBLE_MICROP
           nce = (nc(i,k) + npccn2(i,k)*deltat)
#else
           nce=nc(i,k)
#endif
           nie=(ni(i,k)+nnuccd(i,k)*deltat)

           ! conservation of qc
           !-------------------------------------------------------------------

           dum = (prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ &
                psacws(i,k)+bergs(i,k))*lcldm(i,k)*deltat

           if (dum.gt.qce) then
              ratio = qce/deltat/lcldm(i,k)/(prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+psacws(i,k)+bergs(i,k))*omsm 

              prc(i,k) = prc(i,k)*ratio
              pra(i,k) = pra(i,k)*ratio
              mnuccc(i,k) = mnuccc(i,k)*ratio
              mnucct(i,k) = mnucct(i,k)*ratio  
              msacwi(i,k) = msacwi(i,k)*ratio  
              psacws(i,k) = psacws(i,k)*ratio
              bergs(i,k) = bergs(i,k)*ratio
           end if

           ! conservation of nc
           !-------------------------------------------------------------------
#ifdef GFDL_COMPATIBLE_MICROP
         if (.not. clubb_active) then
           dum = (nprc1(i,k) + npra(i,k) + nnuccc(i,k) + nnucct(i,k) + &
                npsacws(i,k) - nsubc(i,k) - nerosc(i,k))*lcldm(i,k)*deltat
         else
#endif
           dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ &
                npsacws(i,k)-nsubc(i,k))*lcldm(i,k)*deltat
#ifdef GFDL_COMPATIBLE_MICROP
         endif
#endif

           if (dum.gt.nce) then
#ifdef GFDL_COMPATIBLE_MICROP
             if (.not. clubb_active) then
              ratio = nce/deltat/   &
                 ((nprc1(i,k) + npra(i,k) + nnuccc(i,k) + nnucct(i,k) +  &
                 npsacws(i,k) - nsubc(i,k) - nerosc(i,k))*lcldm(i,k))*omsm
             else
#endif
              ratio = nce/deltat/   &
                 ((nprc1(i,k) + npra(i,k) + nnuccc(i,k) + nnucct(i,k) +  &
                 npsacws(i,k) - nsubc(i,k) )*lcldm(i,k))*omsm
#ifdef GFDL_COMPATIBLE_MICROP
             endif
#endif

              nprc1(i,k) = nprc1(i,k)*ratio
              npra(i,k) = npra(i,k)*ratio
              nnuccc(i,k) = nnuccc(i,k)*ratio
              nnucct(i,k) = nnucct(i,k)*ratio
              npsacws(i,k) = npsacws(i,k)*ratio
              nsubc(i,k)=nsubc(i,k)*ratio
#ifdef GFDL_COMPATIBLE_MICROP
              if( .not. clubb_active )  &
              nerosc(i,k)=nerosc(i,k)*ratio
#endif
           end if

           if (do_cldice) then

              ! conservation of qi
              !-------------------------------------------------------------------
              dum = ((-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k))*lcldm(i,k)+(prci(i,k)+ &
                   prai(i,k))*icldm(i,k))*deltat

              if (dum.gt.qie) then

                 ratio = (qie/deltat+(mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k))/ &
                      ((prci(i,k)+prai(i,k))*icldm(i,k))*omsm
                 prci(i,k) = prci(i,k)*ratio
                 prai(i,k) = prai(i,k)*ratio
              end if

              ! conservation of ni
              !-------------------------------------------------------------------
#ifdef GFDL_COMPATIBLE_MICROP
          if( .not. clubb_active ) then
           dum = ((-nnucct(i,k) - nsacwi(i,k))*lcldm(i,k) + (nprci(i,k) + &
                  nprai(i,k) - nsubi(i,k) - nerosi(i,k))*icldm(i,k))*deltat
          else
#endif
              dum = ((-nnucct(i,k)-nsacwi(i,k))*lcldm(i,k)+(nprci(i,k)+ &
                   nprai(i,k)-nsubi(i,k))*icldm(i,k))*deltat
#ifdef GFDL_COMPATIBLE_MICROP
          endif
#endif

              if (dum.gt.nie) then

#ifdef GFDL_COMPATIBLE_MICROP
        if (.not. clubb_active) then
           ratio = (nie/deltat + (nnucct(i,k) + nsacwi(i,k))*lcldm(i,k))/ &
                   ((nprci(i,k) + nprai(i,k) - nsubi(i,k) - nerosi(i,k))* &
                                                          icldm(i,k))*omsm
        else
#endif
            ratio = (nie/deltat+(nnucct(i,k)+nsacwi(i,k))*lcldm(i,k))/ &  
                  ((nprci(i,k)+nprai(i,k)-nsubi(i,k))*icldm(i,k))*omsm
#ifdef GFDL_COMPATIBLE_MICROP
        endif
#endif

                 nprci(i,k) = nprci(i,k)*ratio
                 nprai(i,k) = nprai(i,k)*ratio
                 nsubi(i,k) = nsubi(i,k)*ratio
#ifdef GFDL_COMPATIBLE_MICROP
             if( .not. clubb_active )  &
              nerosi(i,k) = nerosi(i,k)*ratio
#endif
              end if
           end if

           ! for precipitation conservation, use logic that vertical integral 
           ! of tendency from current level to top of model (i.e., qrtot) cannot be negative

           ! conservation of rain mixing rat
           !-------------------------------------------------------------------
           if (((prc(i,k)+pra(i,k))*lcldm(i,k)+(-mnuccr(i,k)+pre(i,k)-pracs(i,k))*&
                cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot(i).lt.0._r8) then

              if (-pre(i,k)+pracs(i,k)+mnuccr(i,k).ge.qsmall) then

                 ratio = (qrtot(i)/(dz(i,k)*rho(i,k))+(prc(i,k)+pra(i,k))*lcldm(i,k))/&
                      ((-pre(i,k)+pracs(i,k)+mnuccr(i,k))*cldmax(i,k))*omsm 

                 pre(i,k) = pre(i,k)*ratio
                 pracs(i,k) = pracs(i,k)*ratio
                 mnuccr(i,k) = mnuccr(i,k)*ratio
              end if
           end if

           ! conservation of nr
           !-------------------------------------------------------------------
           ! for now neglect evaporation of nr
           nsubr(i,k)=0._r8

!--> cjg: modifications incorporated from Huan's code
           if (allow_rain_num_evap)  then
              if (qric(i,k) .ge. qsmall) nsubr(i,k)= max(pre(i,k)/qric(i,k)*  &
                         nric(i,k), -nric(i,k)/deltat)
             endif
!<--cjg
           if ((nprc(i,k)*lcldm(i,k)+(-nnuccr(i,k)+nsubr(i,k)-npracs(i,k)&
                +nragg(i,k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot(i).lt.0._r8) then

              if (-nsubr(i,k)-nragg(i,k)+npracs(i,k)+nnuccr(i,k).ge.qsmall) then

                 ratio = (nrtot(i)/(dz(i,k)*rho(i,k))+nprc(i,k)*lcldm(i,k))/&
                      ((-nsubr(i,k)-nragg(i,k)+npracs(i,k)+nnuccr(i,k))*cldmax(i,k))*omsm

                 nsubr(i,k) = nsubr(i,k)*ratio
                 npracs(i,k) = npracs(i,k)*ratio
                 nnuccr(i,k) = nnuccr(i,k)*ratio
                 nragg(i,k) = nragg(i,k)*ratio
              end if
           end if

           ! conservation of snow mix ratio
           !-------------------------------------------------------------------
           if (((bergs(i,k)+psacws(i,k))*lcldm(i,k)+(prai(i,k)+prci(i,k))*icldm(i,k)+(pracs(i,k)+&
                mnuccr(i,k)+prds(i,k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot(i).lt.0._r8) then

              if (-prds(i,k).ge.qsmall) then

                 ratio = (qstot(i)/(dz(i,k)*rho(i,k))+(bergs(i,k)+psacws(i,k))*lcldm(i,k)+(prai(i,k)+prci(i,k))*icldm(i,k)+&
                      (pracs(i,k)+mnuccr(i,k))*cldmax(i,k))/(-prds(i,k)*cldmax(i,k))*omsm

                 prds(i,k) = prds(i,k)*ratio
              else
                 prds(i,k) = 0._r8
              end if
           end if

           ! conservation of ns
           !-------------------------------------------------------------------
           ! calculate loss of number due to sublimation
           ! for now neglect sublimation of ns
           nsubs(i,k)=0._r8

           if ((nprci(i,k)*icldm(i,k)+(nnuccr(i,k)+nsubs(i,k)+nsagg(i,k))*cldmax(i,k))*&
                dz(i,k)*rho(i,k)+nstot(i).lt.0._r8) then

              if (-nsubs(i,k)-nsagg(i,k).ge.qsmall) then

                 ratio = (nstot(i)/(dz(i,k)*rho(i,k))+nprci(i,k)*icldm(i,k)+&
                      nnuccr(i,k)*cldmax(i,k))/((-nsubs(i,k)-nsagg(i,k))*cldmax(i,k))*omsm

                 nsubs(i,k) = nsubs(i,k)*ratio
                 nsagg(i,k) = nsagg(i,k)*ratio
              end if
           end if

           ! get tendencies due to microphysical conversion processes
           !==========================================================
           ! note: tendencies are multiplied by appropriate cloud/precip 
           ! fraction to get grid-scale values
           ! note: cmei is already grid-average values
#ifdef GFDL_COMPATIBLE_MICROP
          if (.not. clubb_active) then
            qvlat(i,k) = qvlat(i,k) - &
                          (pre(i,k) + prds(i,k))*cldmax(i,k) - cmel(i,k) -  &
                           cmei(i,k) - D_eros_l(i,k) - D_eros_i(i,k)
            tlat(i,k) = tlat(i,k) + ((pre(i,k)*cldmax(i,k) + cmel(i,k) +  &
                          D_eros_l(i,k))*xxlv + (prds(i,k)*cldmax(i,k) +  &
                          cmei(i,k) + D_eros_i(i,k))*xxls + &
                          ((bergs(i,k) + psacws(i,k) + mnuccc(i,k) +   &
                             mnucct(i,k) + msacwi(i,k))*lcldm(i,k) +    &
                             (mnuccr(i,k) + pracs(i,k))*cldmax(i,k) +   &
                                                            berg(i,k))*xlf)
            qctend(i,k) = qctend(i,k) + &
                            (-pra(i,k) - prc(i,k) - mnuccc(i,k) - mnucct(i,k) -   &
                              msacwi(i,k) - psacws(i,k) - bergs(i,k))*   &
                                                           lcldm(i,k) +   &
                                    cmel(i,k) - berg(i,k) + D_eros_l(i,k)
           if (do_cldice) then
            qitend(i,k) = qitend(i,k) + &
                            (mnuccc(i,k) + mnucct(i,k) + msacwi(i,k))*  &
                                                           lcldm(i,k) +   &
                             (-prci(i,k) - prai(i,k))*icldm(i,k) +   &
                                     cmei(i,k) + berg(i,k) + D_eros_i(i,k)
           endif
       else
#endif

           qvlat(i,k) = qvlat(i,k)-(pre(i,k)+prds(i,k))*cldmax(i,k)-cmei(i,k) 

           tlat(i,k) = tlat(i,k)+((pre(i,k)*cldmax(i,k)) &
                *xxlv+(prds(i,k)*cldmax(i,k)+cmei(i,k))*xxls+ &
                ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(mnuccr(i,k)+ &
                pracs(i,k))*cldmax(i,k)+berg(i,k))*xlf)

           qctend(i,k) = qctend(i,k)+ &
                (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- & 
                psacws(i,k)-bergs(i,k))*lcldm(i,k)-berg(i,k)

           if (do_cldice) then
              qitend(i,k) = qitend(i,k)+ &
                   (mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(-prci(i,k)- & 
                   prai(i,k))*icldm(i,k)+cmei(i,k)+berg(i,k)
           end if
#ifdef GFDL_COMPATIBLE_MICROP
       end if
#endif

           qrtend(i,k) = qrtend(i,k)+ &
                (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
                mnuccr(i,k))*cldmax(i,k)

           qstend(i,k) = qstend(i,k)+ &
                (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(prds(i,k)+ &
                pracs(i,k)+mnuccr(i,k))*cldmax(i,k)

           ! add output for cmei (accumulate)
           cmeitot(i,k) = cmeitot(i,k) + cmei(i,k)

           ! assign variables for trop_mozart, these are grid-average
           !-------------------------------------------------------------------
           ! evaporation/sublimation is stored here as positive term

           evapsnow(i,k) = evapsnow(i,k)-prds(i,k)*cldmax(i,k)
           nevapr(i,k) = nevapr(i,k)-pre(i,k)*cldmax(i,k)

           ! change to make sure prain is positive: do not remove snow from
           ! prain used for wet deposition
           prain(i,k) = prain(i,k)+(pra(i,k)+prc(i,k))*lcldm(i,k)+(-pracs(i,k)- &
                mnuccr(i,k))*cldmax(i,k)
           prodsnow(i,k) = prodsnow(i,k)+(prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(&
                pracs(i,k)+mnuccr(i,k))*cldmax(i,k)

           ! following are used to calculate 1st order conversion rate of cloud water
           !    to rain and snow (1/s), for later use in aerosol wet removal routine
           ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
           !    used to calculate pra, prc, ... in this routine
           ! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
           !                      (no cloud ice or bergeron terms)
           ! qcsum_rate1ord     = sum over iterations{ qc used in calculation of the transfer terms }

           qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) + (pra(i,k)+prc(i,k)+psacws(i,k))*lcldm(i,k) 
           qcsum_rate1ord(i,k) = qcsum_rate1ord(i,k) + qc(i,k) 

           ! microphysics output, note this is grid-averaged
#ifdef GFDL_COMPATIBLE_MICROP
            preo(i,k) = preo(i,k) + pre(i,k)*cldmax(i,k)
            prdso(i,k) = prdso(i,k) + prds(i,k)*cldmax(i,k)
          if (.not. clubb_active) then
            cmelo(i,k) = cmelo(i,k) + cmel(i,k)
            eroslo(i,k) = eroslo(i,k) + D_eros_l(i,k)
            erosio(i,k) = erosio(i,k) + D_eros_i(i,k)
          endif
#endif
           pratot(i,k)=pratot(i,k)+pra(i,k)*lcldm(i,k)
           prctot(i,k)=prctot(i,k)+prc(i,k)*lcldm(i,k)
           mnuccctot(i,k)=mnuccctot(i,k)+mnuccc(i,k)*lcldm(i,k)
           mnuccttot(i,k)=mnuccttot(i,k)+mnucct(i,k)*lcldm(i,k)
           mnuccdtot(i,k)=mnuccdtot(i,k)+mnuccd(i,k)*lcldm(i,k)
           msacwitot(i,k)=msacwitot(i,k)+msacwi(i,k)*lcldm(i,k)
           psacwstot(i,k)=psacwstot(i,k)+psacws(i,k)*lcldm(i,k)
           bergstot(i,k)=bergstot(i,k)+bergs(i,k)*lcldm(i,k)
           bergtot(i,k)=bergtot(i,k)+berg(i,k)
           prcitot(i,k)=prcitot(i,k)+prci(i,k)*icldm(i,k)
           praitot(i,k)=praitot(i,k)+prai(i,k)*icldm(i,k)
           mnuccrtot(i,k)=mnuccrtot(i,k)+mnuccr(i,k)*cldmax(i,k)
           pracstot(i,k)=pracstot(i,k)+pracs(i,k)*cldmax(i,k)
#ifdef GFDL_COMPATIBLE_MICROP
         if (.not. clubb_active) then
           nctend(i,k) = nctend(i,k)+ npccn2(i,k)        + &
                (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) & 
                + nerosc(i,k) -npra(i,k)-nprc1(i,k))*lcldm(i,k)

           if (do_cldice) then
              nitend(i,k) = nitend(i,k)+ nnuccd(i,k)+ &
                   (nnucct(i,k)+nsacwi(i,k))*lcldm(i,k)+(nsubi(i,k) + &
                      nerosi(i,k) -nprci(i,k)- &
                   nprai(i,k))*icldm(i,k)
           end if
        else
#endif
           nctend(i,k) = nctend(i,k)+&
                (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) & 
                -npra(i,k)-nprc1(i,k))*lcldm(i,k)

           if (do_cldice) then
              nitend(i,k) = nitend(i,k)+ nnuccd(i,k)+ &
                   (nnucct(i,k)+nsacwi(i,k))*lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- &
                   nprai(i,k))*icldm(i,k)
            endif
#ifdef GFDL_COMPATIBLE_MICROP
        endif
#endif

           nstend(i,k) = nstend(i,k)+(nsubs(i,k)+ &
                nsagg(i,k)+nnuccr(i,k))*cldmax(i,k)+nprci(i,k)*icldm(i,k)

           nrtend(i,k) = nrtend(i,k)+ &
                nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
                +nragg(i,k))*cldmax(i,k)

#ifdef GFDL_COMPATIBLE_MICROP
! save current tendencies so that upcoming adjustment may be captured
            IF (diag_id%qndt_nucclim +     &
                                      diag_id%qn_nucclim_col  > 0) THEN
              nucclim(i,k) = nctend(i,k)
             END IF
 
            IF (diag_id%qnidt_nucclim1 +    &
                                       diag_id%qni_nucclim1_col > 0) THEN
              nucclim1i(i,k) = nitend(i,k)
            END IF
#endif
           ! make sure that ni at advanced time step does not exceed
           ! maximum (existing N + source terms*dt), which is possible if mtime < deltat
           ! note that currently mtime = deltat
           !================================================================
#ifdef GFDL_COMPATIBLE_MICROP
          if (nctend(i,k) .gt. 0._r8 .and.    &
                   nc(i,k) + nctend(i,k)*deltat .gt. ncmax(i,k)) then
              nctend(i,k) = max(0._r8, (ncmax(i,k) - nc(i,k))/deltat)
            end if
#endif

           if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax(i,k)) then
              nitend(i,k)=max(0._r8,(nimax(i,k)-ni(i,k))/deltat)
           end if

#ifdef GFDL_COMPATIBLE_MICROP
! complete diagnostic calculation 
            IF (diag_id%qndt_nucclim +     &
                                     diag_id%qn_nucclim_col  > 0) THEN
              nucclim(i,k) = nctend(i,k) - nucclim(i,k)
            END IF
            IF (diag_id%qnidt_nucclim1 +     &
                                    diag_id%qni_nucclim1_col > 0) THEN
              nucclim1i(i,k) = nitend(i,k) - nucclim1i(i,k)
            END IF
#endif
        end do
           
        ! End of "administration" loop

        ! get final values for precipitation q and N, based on
        ! flux of precip from above, source/sink term, and terminal fallspeed
        ! see eq. 15-16 in MG2008

        do i = 1, mgncol

           ! rain

           if (qric(i,k).ge.qsmall) then
              if (k .eq. 1) then
                 qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(i,k)
                 nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(i,k)
              else
                 qric(i,k) = (rho(i,k-1)*umr(i,k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
                      (rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(i,k)*rho(i,k)*cldmax(i,k))
                 nric(i,k) = (rho(i,k-1)*unr(i,k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
                      (rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(i,k)*rho(i,k)*cldmax(i,k))

              end if
           else
              qric(i,k)=0._r8
              nric(i,k)=0._r8
           end if

           ! snow

           if (qsic(i,k).ge.qsmall) then
              if (k .eq. 1) then
                 qsic(i,k)=qstend(i,k)*dz(i,k)/cldmax(i,k)/ums(i,k)
                 nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(i,k)
              else
                 qsic(i,k) = (rho(i,k-1)*ums(i,k-1)*qsic(i,k-1)*cldmax(i,k-1)+ &
                      (rho(i,k)*dz(i,k)*qstend(i,k)))/(ums(i,k)*rho(i,k)*cldmax(i,k))
                 nsic(i,k) = (rho(i,k-1)*uns(i,k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
                      (rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(i,k)*rho(i,k)*cldmax(i,k))
              end if
           else
              qsic(i,k)=0._r8
              nsic(i,k)=0._r8
           end if

           ! calculate precipitation flux at surface
           !=========================================================
           ! divide by density of water to get units of m/s

           prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+&
                qstend(i,k)*dz(i,k)*rho(i,k))/rhow
           preci(i) = preci(i)+qstend(i,k)*dz(i,k)*rho(i,k)/rhow

           ! convert rain rate from m/s to mm/hr

           rainrt(i,k)=rainrt(i,k) + (qric(i,k)*rho(i,k)*umr(i,k)/rhow*3600._r8*1000._r8)

           ! vertically-integrated precip source/sink terms (note: grid-averaged)

           qrtot(i) = max(qrtot(i)+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
           qstot(i) = max(qstot(i)+qstend(i,k)*dz(i,k)*rho(i,k),0._r8)
           nrtot(i) = max(nrtot(i)+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
           nstot(i) = max(nstot(i)+nstend(i,k)*dz(i,k)*rho(i,k),0._r8)

           ! calculate melting and freezing of precip
           !=========================================================

           ! melt snow at +2 C

           ! Difference in amount of heat between temperature at which
           ! all snow melts, and the current state.
           dum = (snowmelt - t(i,k) - tlat(i,k)/cpp*deltat) * cpp

           ! Test if temperature is above threshold and snow is present.
           if (dum < 0._r8 .and. qstot(i) > 0._r8) then

              ! (negative) heat produced if all snow is melted
              dum1 = -xlf * qstot(i)/(dz(i,k)*rho(i,k))*deltat

              ! ratio of heating needed to get to the threshold, to
              ! the total heating from melting everything, is equal
              ! to the proportion of snow that actually melts
              dum = min(1._r8, dum/dum1)
              dum = max(0._r8, dum)

              ! Melt snow
              qric(i,k)=qric(i,k)+dum*qsic(i,k)
              nric(i,k)=nric(i,k)+dum*nsic(i,k)
              qsic(i,k)=(1._r8-dum)*qsic(i,k)
              nsic(i,k)=(1._r8-dum)*nsic(i,k)

#ifdef GFDL_COMPATIBLE_MICROP
                if (diag_id%snow_melt + diag_id%snow_melt_col > 0) &
                  diag_4l(i,j,k, diag_pt%snow_melt) =   &
                         diag_4l(i,j,k, diag_pt%snow_melt) +   &
                                 dum*preci(i)*rhow/(rho(i,k)*dz(i,k))
#endif

              qrtot(i)=qrtot(i)+dum*qstot(i)
              nrtot(i)=nrtot(i)+dum*nstot(i)
              qstot(i)=(1._r8-dum)*qstot(i)
              nstot(i)=(1._r8-dum)*nstot(i)

              preci(i)=(1._r8-dum)*preci(i)

              ! Get heating tendency based on proportion of snow that
              ! actually melts.
              dum1 = dum * dum1/deltat

              meltsdttot(i,k)=meltsdttot(i,k) + dum1
              tlat(i,k)=tlat(i,k)+dum1

           end if

           ! freeze all rain at -5C for Arctic
           !=========================================================

           ! Difference in amount of heat between temperature at which
           ! all rain freezes, and the current state.
           dum = (rainfrze - t(i,k) - tlat(i,k)/cpp*deltat) * cpp

           ! Test if temperature is below threshold and snow is present.
           if (dum > 0._r8 .and. qrtot(i) > 0._r8) then

              ! heat produced if all rain freezes
              dum1 = xlf * qrtot(i)/(dz(i,k)*rho(i,k))*deltat

              ! ratio of heating needed to get to the threshold, to
              ! the total heating from freezing everything, is equal
              ! to the proportion of rain that actually freezes
              dum = min(1._r8, dum/dum1)
              dum = max(0._r8, dum)

#ifdef GFDL_COMPATIBLE_MICROP
                diag_4l(i,j,k, diag_pt%rain_freeze) =  &
                         diag_4l(i,j,k, diag_pt%rain_freeze) +  &
                                            dum*(prect(i) - preci(i))*  &
                                                   rhow /(rho(i,k)*dz(i,k))
#endif
              ! Freeze rain
              qsic(i,k)=qsic(i,k)+dum*qric(i,k)
              nsic(i,k)=nsic(i,k)+dum*nric(i,k)
              qric(i,k)=(1._r8-dum)*qric(i,k)
              nric(i,k)=(1._r8-dum)*nric(i,k)

              qstot(i)=qstot(i)+dum*qrtot(i)
              nstot(i)=nstot(i)+dum*nrtot(i)
              qrtot(i)=(1._r8-dum)*qrtot(i)
              nrtot(i)=(1._r8-dum)*nrtot(i)

              preci(i)=preci(i)+dum*(prect(i)-preci(i))

              ! Get heating tendency based on proportion of rain that
              ! actually freezes.
              dum1 = dum * dum1/deltat

              frzrdttot(i,k)=frzrdttot(i,k) + dum1
              tlat(i,k)=tlat(i,k)+dum1

           end if

           ! if rain/snow mix ratio is zero so should number concentration
           !=========================================================
           
           if (qsic(i,k) < qsmall) then
              qsic(i,k)=0._r8
              nsic(i,k)=0._r8
           end if

           if (qric(i,k) < qsmall) then
              qric(i,k)=0._r8
              nric(i,k)=0._r8
           end if

           ! make sure number concentration is a positive number to avoid 
           ! taking root of negative

           nric(i,k)=max(nric(i,k),0._r8)
           nsic(i,k)=max(nsic(i,k),0._r8)

           ! get size distribution parameters for fallspeed calculations
           !=========================================================

           ! rain

           call size_dist_param_rain(qric(i,k), nric(i,k), lamr(i,k), n0r(i,k))

           if (lamr(i,k).ge.qsmall) then

              ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)

              unr(i,k) = min(arn(i,k)*cons4/lamr(i,k)**br,9.1_r8*rhof(i,k))
              umr(i,k) = min(arn(i,k)*cons5/(6._r8*lamr(i,k)**br),9.1_r8*rhof(i,k))

           else
              umr(i,k)=0._r8
              unr(i,k)=0._r8
           end if

           !......................................................................
           ! snow

           call size_dist_param_snow(qsic(i,k), nsic(i,k), lams(i,k), n0s(i,k))

           if (lams(i,k) > 0._r8) then

              ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)

              ums(i,k) = min(asn(i,k)*cons8/(6._r8*lams(i,k)**bs),1.2_r8*rhof(i,k))
              uns(i,k) = min(asn(i,k)*cons7/lams(i,k)**bs,1.2_r8*rhof(i,k))

           else
              ums(i,k) = 0._r8
              uns(i,k) = 0._r8
           end if

        end do

        ! Done with vertical dependencies from precipitation.

     end do micro_vert_loop ! end k loop

     ! sum over sub-step for average process rates
     !-----------------------------------------------------
     ! convert rain/snow q and N for output to history, note, 
     ! output is for gridbox average

     ! calculate precip fluxes and adding them to summing sub-stepping variables
     ! calculate the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s)
     ! ---------------------------------------------------------------------

     rflx(:,2:) = rflx(:,2:) + (qric*rho*umr*cldmax)

     dumr(:,1)  = int_to_mid(:,1) * qric(:,1)
     qrout(:,1)  = qrout(:,1)  + int_to_mid(:,1)*qric(:,1)*cldmax(:,1)

     dumr(:,2:)  = interp_to_mid(qric,int_to_mid(:,2:))
     qrout(:,2:) = qrout(:,2:) + interp_to_mid(qric * cldmax, int_to_mid(:,2:))

     dumnr(:,1) = int_to_mid(:,1) * nric(:,1)
     nrout(:,1)  = nrout(:,1)  + &
          (int_to_mid(:,1)*nric(:,1)*cldmax(:,1)*rho(:,1))

     dumnr(:,2:) = interp_to_mid(nric,int_to_mid(:,2:))
     nrout(:,2:) = nrout(:,2:) + interp_to_mid(nric * cldmax * rho, int_to_mid(:,2:))
     ! Calculate rercld

     ! calculate mean size of combined rain and cloud water
 
     ! hm 11-22-11 modify to interpolate rain from interface to mid-point
     ! logic is to interpolate rain mass and number, then recalculate PSD
     ! parameters to get relevant parameters for mean size

     ! interpolate rain mass and number, store in dummy variables
     
     ! calculate n0r and lamr from interpolated mid-point rain mass and number
     ! divide by precip fraction to get in-precip (local) values of
     ! rain mass and number, divide by rhow to get rain number in kg^-1

     call size_dist_param_rain(dumr, dumnr, lamr, n0r)

#ifdef GFDL_COMPATIBLE_MICROP
     do k=1,nlev
       do i=1,mgncol
     call calc_rercld(lamr(i,k), n0r(i,k), lamc(i,k), cdist1(i,k), pgam(i,k), dumr(i,k), qcic(i,k), &
          arcld(i,k), rercld(i,k))
    end do
    end do
#else
     call calc_rercld(lamr, n0r, lamc, cdist1, pgam, dumr, qcic, &
          arcld, rercld)
#endif

     nsout(:,1)  = nsout(:,1)  + &
          (int_to_mid(:,1)*nsic(:,1)*cldmax(:,1)*rho(:,1))
     nsout(:,2:) = nsout(:,2:) + interp_to_mid(nsic * cldmax * rho, int_to_mid(:,2:))

     qsout(:,1)  = qsout(:,1)  + int_to_mid(:,1)*qsic(:,1)*cldmax(:,1)
     qsout(:,2:) = qsout(:,2:) + interp_to_mid(qsic * cldmax, int_to_mid(:,2:))

     sflx(:,2:) = sflx(:,2:) + (qsic*rho*ums*cldmax)

     ! Sum into other variables that accumulate over substeps.
     tlat1 = tlat1 + tlat
     t = t + tlat*deltat/cpp

     qvlat1 = qvlat1 + qvlat
     q = q + qvlat*deltat

     qctend1 = qctend1 + qctend
     qc = qc + qctend*deltat

     qitend1 = qitend1 + qitend
     qi = qi + qitend*deltat

     nctend1 = nctend1 + nctend
     nc = nc + nctend*deltat

     nitend1 = nitend1 + nitend
     ni = ni + nitend*deltat

     prect1 = prect1 + prect
     preci1 = preci1 + preci

#ifdef GFDL_COMPATIBLE_MICROP
!droplet number
!           npccno     = npccno   + npccn*mtime
            npccno     = npccno   + npccn2
            nnuccco    = nnuccco - nnuccc*lcldm
            nnuccto    = nnuccto  - nnucct*lcldm
            npsacwso   = npsacwso - npsacws*lcldm
            nsubco     = nsubco   + nsubc*lcldm
          if (.not. clubb_active) &
             nerosco    = nerosco  + nerosc*lcldm
            nprao      = nprao    - npra*lcldm
             nprc1o     = nprc1o   - nprc1*lcldm
             nucclimo   = nucclimo + nucclim
 
!ice number
!           nnuccdo    = nnuccdo + nnuccd*mtime
            nnuccdo    = nnuccdo + nnuccd
            nsacwio    = nsacwio + nsacwi*lcldm
            nsubio     = nsubio  + nsubi*icldm
          if (.not. clubb_active) &
            nerosio    = nerosio + nerosi*icldm
            nprcio     = nprcio  - nprci*icldm
            npraio     = npraio  - nprai*icldm
            nucclim1io = nucclim1io + nucclim1i
#endif

  end do substepping ! it loop, sub-step

  ! divide rain radius over substeps for average
  where (arcld > 0) rercld = rercld/arcld

  ! convert dt from sub-step back to full time step
  !-------------------------------------------------------------------
  deltat = deltatin

  ! assign variables back to start-of-timestep values before updating after sub-steps 
  !================================================================================
  
  call pack_array(qn, mgcols, top_lev, q)
  call pack_array(tn, mgcols, top_lev, t)
  call pack_array(qcn, mgcols, top_lev, qc)
  call pack_array(ncn, mgcols, top_lev, nc)
  call pack_array(qin, mgcols, top_lev, qi)
  call pack_array(nin, mgcols, top_lev, ni)

  !.............................................................................

  ! divide precip rate by number of sub-steps to get average over time step
  
  prect = prect1/real(iter)
  preci = preci1/real(iter)

#ifdef GFDL_COMPATIBLE_MICROP
! convert unit  from m/s to kg/m2/s by multiply water density 1000 kg/m3
        diag_4l(:,j,:, diag_pt%snow_melt) =  &
                          diag_4l(:,j,:, diag_pt%snow_melt)/real(iter)
        diag_4l(:,j,:, diag_pt%rain_freeze) =  &
                          diag_4l(:,j,:, diag_pt%rain_freeze)/real(iter)
! <--- h1g, 2010-11-15
#endif
  ! divide microphysical tendencies by number of sub-steps to get average over time step
  !================================================================================

  tlat = tlat1/real(iter)
  qvlat = qvlat1/real(iter)
  qctend = qctend1/real(iter)
  qitend = qitend1/real(iter)
  nctend = nctend1/real(iter)
  nitend = nitend1/real(iter)

  ! Re-apply droplet activation tendency
#ifndef GFDL_COMPATIBLE_MICROP
  nctend = nctend + npccn
#endif

  rainrt = rainrt/real(iter)

  ! divide by number of sub-steps to find final values
  rflx = rflx/real(iter)
  sflx = sflx/real(iter)

  ! divide output precip q and N by number of sub-steps to get average over time step
  !================================================================================

  qrout = qrout/real(iter)
  qsout = qsout/real(iter)
  nrout = nrout/real(iter)
  nsout = nsout/real(iter)

  ! divide trop_mozart variables by number of sub-steps to get average over time step 
  !================================================================================

  nevapr = nevapr/real(iter)
  evapsnow = evapsnow/real(iter)
  prain = prain/real(iter)
  prodsnow = prodsnow/real(iter)

  ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
  nevapr = nevapr + evapsnow
  prain = prain + prodsnow

  cmeout = cmeout/real(iter)

  cmeitot = cmeitot/real(iter)
  meltsdttot = meltsdttot/real(iter)
  frzrdttot  = frzrdttot /real(iter)

  ! microphysics output
#ifdef GFDL_COMPATIBLE_MICROP
! vapor
          preo = preo/real(iter)
          prdso = prdso/real(iter)
          if (.not. clubb_active) then
          cmelo =  cmelo/real(iter)
          eroslo = eroslo/real(iter)
          erosio = erosio/real(iter)
          endif
#endif
  pratot=pratot/real(iter)
  prctot=prctot/real(iter)
  mnuccctot=mnuccctot/real(iter)
  mnuccttot=mnuccttot/real(iter)
  msacwitot=msacwitot/real(iter)
  psacwstot=psacwstot/real(iter)
  bergstot=bergstot/real(iter)
  bergtot=bergtot/real(iter)
  prcitot=prcitot/real(iter)
  praitot=praitot/real(iter)

  mnuccrtot=mnuccrtot/real(iter)
  pracstot =pracstot /real(iter)

  mnuccdtot=mnuccdtot/real(iter)

  sed_col_loop: do i=1,mgncol

     do k=1,nlev

        ! calculate sedimentation for cloud water and ice
        !================================================================================

        ! update in-cloud cloud mixing ratio and number concentration 
        ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
        ! note: these are in-cloud values***, hence we divide by cloud fraction

        dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
        dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
        dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
        dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)

        ! hm add 6/2/11 switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k)=ncnst/rho(i,k)
        end if

        ! hm add 6/2/11 switch for specification of cloud ice number
        if (nccons  .or. liq_sedi_mass_only) then       ! h1g, 2016-03-07
           dumni(i,k)=ninst/rho(i,k)
        end if

        ! obtain new slope parameter to avoid possible singularity

        call size_dist_param_ice(dumi(i,k), dumni(i,k), lami(i,k), n0i(i,k))

        call size_dist_param_liq(dumc(i,k), dumnc(i,k), cdnl, rho(i,k), .true., &
             pgam(i,k), lamc(i,k))

        ! calculate number and mass weighted fall velocity for droplets and cloud ice
        !-------------------------------------------------------------------


        if (dumc(i,k).ge.qsmall) then

           vtrmc(i,k)=acn(i,k)*gamma(4._r8+bc+pgam(i,k))/ &
                (lamc(i,k)**bc*gamma(pgam(i,k)+4._r8))

           fc(k) = g*rho(i,k)*vtrmc(i,k)

           fnc(k) = g*rho(i,k)* &
                acn(i,k)*gamma(1._r8+bc+pgam(i,k))/ &
                (lamc(i,k)**bc*gamma(pgam(i,k)+1._r8))

 !---> h1g, 2016-03-07
           if ( no_liq_sedi ) then
              fc(k)  = 0.0
              fnc(k) = 0.0
           endif
!<--- h1g, 2016-03-07

        else
           fc(k) = 0._r8
           fnc(k)= 0._r8
        end if

        ! calculate number and mass weighted fall velocity for cloud ice

        if (dumi(i,k).ge.qsmall) then

           vtrmi(i,k)=min(ain(i,k)*cons17/(6._r8*lami(i,k)**bi), &
                1.2_r8*rhof(i,k))

           fi(k) = g*rho(i,k)*vtrmi(i,k)
           fni(k) = g*rho(i,k)* &
                min(ain(i,k)*cons16/lami(i,k)**bi,1.2_r8*rhof(i,k))
#ifdef GFDL_COMPATIBLE_MICROP
          if (diag_id%vfall > 0) diag_4l(i,j,k,diag_pt%vfall) = vtrmi(i,k)
#endif
        else
           fi(k) = 0._r8
           fni(k)= 0._r8
#ifdef GFDL_COMPATIBLE_MICROP
          if (diag_id%vfall > 0) diag_4l(i,j,k,diag_pt%vfall) = 0._r8
#endif
        end if

        ! redefine dummy variables - sedimentation is calculated over grid-scale
        ! quantities to ensure conservation

        dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
        dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
        dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)

        if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8

     end do       !!! vertical loop

     ! initialize nstep for sedimentation sub-steps

     ! calculate number of split time steps to ensure courant stability criteria
     ! for sedimentation calculations
     !-------------------------------------------------------------------
     nstep = 1 + int(max( &
          maxval( fi/pdel(i,:)), &
          maxval( fc/pdel(i,:)), &
          maxval(fni/pdel(i,:)), &
          maxval(fnc/pdel(i,:))) &
          * deltat)

#ifdef GFDL_COMPATIBLE_MICROP
     diag_4l(i,j,:,diag_pt%qnidt_sedi) = 0.   
     diag_4l(i,j,:,diag_pt%qndt_sedi) = 0.   
#endif

     ! loop over sedimentation sub-time step to ensure stability
     !==============================================================
     do n = 1,nstep

        if (do_cldice) then
           falouti  = fi  * dumi(i,:)
           faloutni = fni * dumni(i,:)
        else
           falouti  = 0._r8
           faloutni = 0._r8
        end if
        faloutc  = fc  * dumc(i,:)
        faloutnc = fnc * dumnc(i,:)

        ! top of model

        k = 1
        faltndi = falouti(k)/pdel(i,k)
        faltndni = faloutni(k)/pdel(i,k)
        faltndc = faloutc(k)/pdel(i,k)
        faltndnc = faloutnc(k)/pdel(i,k)

        ! add fallout terms to microphysical tendencies

        qitend(i,k) = qitend(i,k)-faltndi/nstep
        nitend(i,k) = nitend(i,k)-faltndni/nstep
        qctend(i,k) = qctend(i,k)-faltndc/nstep
        nctend(i,k) = nctend(i,k)-faltndnc/nstep

        ! sedimentation tendencies for output
        qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
        qisedten(i,k)=qisedten(i,k)-faltndi/nstep

#ifdef GFDL_COMPATIBLE_MICROP
          IF ( diag_id%qndt_sedi + diag_id%qn_sedi_col > 0 ) &
                 diag_4l(i,j,k,diag_pt%qndt_sedi) =     &
                       diag_4l(i,j,k,diag_pt%qndt_sedi) - faltndnc/nstep

          IF ( diag_id%qnidt_sedi +  diag_id%qni_sedi_col > 0 ) &
                 diag_4l(i,j,k,diag_pt%qnidt_sedi) =    &
                      diag_4l(i,j,k,diag_pt%qnidt_sedi) - faltndni/nstep
#endif

        dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
        dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
        dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
        dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep

        do k = 2,nlev

           ! for cloud liquid and ice, if cloud fraction increases with height
           ! then add flux from above to both vapor and cloud water of current level
           ! this means that flux entering clear portion of cell from above evaporates
           ! instantly

           dum=lcldm(i,k)/lcldm(i,k-1)
           dum=min(dum,1._r8)
           dum1=icldm(i,k)/icldm(i,k-1)
           dum1=min(dum1,1._r8)

           faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k)
           faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k)
           faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k)
           faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k)
           faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k)
           faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k)

           ! add fallout terms to eulerian tendencies

           qitend(i,k) = qitend(i,k)-faltndi/nstep
           nitend(i,k) = nitend(i,k)-faltndni/nstep
           qctend(i,k) = qctend(i,k)-faltndc/nstep
           nctend(i,k) = nctend(i,k)-faltndnc/nstep

           ! sedimentation tendencies for output
           qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
           qisedten(i,k)=qisedten(i,k)-faltndi/nstep

           ! add terms to to evap/sub of cloud water

#ifdef GFDL_COMPATIBLE_MICROP
            IF ( diag_id%qndt_sedi  + diag_id%qn_sedi_col > 0 ) &
               diag_4l(i,j,k,diag_pt%qndt_sedi) =    &
                   diag_4l(i,j,k,diag_pt%qndt_sedi) - faltndnc/nstep

            IF ( diag_id%qnidt_sedi +  diag_id%qni_sedi_col > 0 ) &
                  diag_4l(i,j,k,diag_pt%qnidt_sedi) =    &
                        diag_4l(i,j,k,diag_pt%qnidt_sedi) - faltndni/nstep
#endif
           qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep
           ! for output
           qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep
           qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep
           ! for output
           qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep

           tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep
           tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep

           dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
           dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
           dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
           dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep

! ---> h1g, 2014-05-23 keep the sedimentation bugs in previous MG versions
           if( .not. SedBugFix ) then   !sedimentation increases with press ure (or decreases with altitude)     
             Fni(K) = MAX(Fni(K)/pdel(i,K), Fni(K-1)/pdel(i,K-1))*pdel(i,K)
             FI(K) = MAX(FI(K)/pdel(i,K), FI(K-1)/pdel(i,K-1))*pdel(i,K)
             fnc(k) = max(fnc(k)/pdel(i,k), fnc(k-1)/pdel(i,k-1))*pdel(i,k)
             Fc(K) = MAX(Fc(K)/pdel(i,K), Fc(K-1)/pdel(i,K-1))*pdel(i,K)
           endif  ! SedBugFix 
! <--- h1g, 2014-05-23

        end do   !! k loop

        ! units below are m/s
        ! cloud water/ice sedimentation flux at surface 
        ! is added to precip flux at surface to get total precip (cloud + precip water)
        ! rate

        prect(i) = prect(i)+(faloutc(nlev)+falouti(nlev))/g/nstep/1000._r8
        preci(i) = preci(i)+(              falouti(nlev))/g/nstep/1000._r8

#ifdef GFDL_COMPATIBLE_MICROP
!  for the diagnostic of surface precip
           k = 1
           IF ( diag_id%sedi_sfc > 0 ) &
                   diag_4l(i,j,1,diag_pt%sedi_sfc) =    &
                          diag_4l(i,j,1,diag_pt%sedi_sfc) +  &
                                            faloutc(nlev)/g/nstep/rhow
          IF ( diag_id%sedi_ice > 0 ) &
                  diag_4l(i,j,1,diag_pt%sedi_ice) =   &
                           diag_4l(i,j,1,diag_pt%sedi_ice) +   &
                                             falouti(nlev)/g/nstep/rhoi
#endif

     end do   !! nstep loop

     ! end sedimentation
     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

     ! get new update for variables that includes sedimentation tendency
     ! note : here dum variables are grid-average, NOT in-cloud

     do k=1,nlev

        dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
        dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)

        ! hm add 6/2/11 switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k)=ncnst/rho(i,k)*lcldm(i,k)
        end if

        ! hm add 6/2/11 switch for specification of cloud ice number
        if (nicons) then
           dumni(i,k)=ninst/rho(i,k)*icldm(i,k)
        end if

        if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8

        ! calculate instantaneous processes (melting, homogeneous freezing)
        !====================================================================

        if (do_cldice) then
           if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then
              if (dumi(i,k) > 0._r8) then

                 ! limit so that melting does not push temperature below freezing
                 !-----------------------------------------------------------------
                 dum = -dumi(i,k)*xlf/cpp
                 if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then
                    dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf
                    dum = dum/dumi(i,k)
                    dum = max(0._r8,dum)
                    dum = min(1._r8,dum)
                 else
                    dum = 1._r8
                 end if

                 qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat

                 ! for output
                 melttot(i,k)=dum*dumi(i,k)/deltat


                 ! assume melting ice produces droplet
                 ! mean volume radius of 8 micron

#ifdef GFDL_COMPATIBLE_MICROP
              IF (diag_id%qndt_melt + diag_id%qn_melt_col > 0) &
                   diag_4l(i,j,k,diag_pt%qndt_melt) =  nctend(i,k)
              IF (diag_id%qidt_melt2  + diag_id%qi_melt2_col  > 0) &
                   diag_4l(i,j,k,diag_pt%qidt_melt2) =  qitend(i,k)
              IF (diag_id%qnidt_melt +  diag_id%qni_melt_col > 0) &
                   diag_4l(i,j,k,diag_pt%qnidt_melt) =  nitend(i,k)
#endif
                 nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ &
                      (4._r8*pi*5.12e-16_r8*rhow)

                 qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat
                 nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat
                 tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat
#ifdef GFDL_COMPATIBLE_MICROP
              IF (diag_id%qndt_melt + diag_id%qn_melt_col > 0) &
                   diag_4l(i,j,k,diag_pt%qndt_melt) =    &
                          nctend(i,k) - diag_4l(i,j,k,diag_pt%qndt_melt)
              IF (diag_id%qidt_melt2  + diag_id%qi_melt2_col  > 0) &
                   diag_4l(i,j,k,diag_pt%qidt_melt2) =    &
                          qitend(i,k) - diag_4l(i,j,k,diag_pt%qidt_melt2)
              IF (diag_id%qnidt_melt +  diag_id%qni_melt_col > 0) &
                   diag_4l(i,j,k,diag_pt%qnidt_melt) =    &
                          nitend(i,k) - diag_4l(i,j,k,diag_pt%qnidt_melt)
#endif
              end if
           end if

           ! homogeneously freeze droplets at -40 C
           !-----------------------------------------------------------------

#ifdef GFDL_COMPATIBLE_MICROP
           if (t(i,k)+tlat(i,k)/cpp*deltat < tmelt -40._r8) then
#else
           if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then
#endif
              if (dumc(i,k) > 0._r8) then

                 ! limit so that freezing does not push temperature above threshold
                 dum = dumc(i,k)*xlf/cpp
#ifdef GFDL_COMPATIBLE_MICROP
                 if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.tmelt-40._r8) then
                    dum = -(t(i,k)+tlat(i,k)/cpp*deltat-(tmelt-40._r8))*cpp/xlf
#else
                 if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then
                    dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf
#endif
                    dum = dum/dumc(i,k)
                    dum = max(0._r8,dum)
                    dum = min(1._r8,dum)
                 else
                    dum = 1._r8
                 end if

                 qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat
                 ! for output
                 homotot(i,k)=dum*dumc(i,k)/deltat

                 ! assume 25 micron mean volume radius of homogeneously frozen droplets
                 ! consistent with size of detrained ice in stratiform.F90

#ifdef GFDL_COMPATIBLE_MICROP
              IF (diag_id%qldt_freez + diag_id%ql_freez_col > 0) &
                    diag_4l(i,j,k,diag_pt%qldt_freez) = qctend(i,k)
               sum_freeze(i,k) = qctend(i,k)
               IF ( diag_id%qndt_ihom  + diag_id%qn_ihom_col > 0) &
                    diag_4l(i,j,k,diag_pt%qndt_ihom) =  nctend(i,k)
#endif

                 nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* &
                      500._r8)/deltat
                 qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat
                 if ( include_homo_freez_in_qn ) &        ! h1g, 2014-07-22
                 nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat
                 tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat

#ifdef GFDL_COMPATIBLE_MICROP
               IF (diag_id%qldt_freez + diag_id%ql_freez_col > 0) &
                    diag_4l(i,j,k,diag_pt%qldt_freez) =    &
                          qctend(i,k) - diag_4l(i,j,k,diag_pt%qldt_freez)
               sum_freeze(i,k) = -(qctend(i,k) - sum_freeze(i,k))
               IF (diag_id%qndt_ihom  + diag_id%qn_ihom_col > 0) &
                    diag_4l(i,j,k,diag_pt%qndt_ihom) =    &
                            nctend(i,k) - diag_4l(i,j,k,diag_pt%qndt_ihom)
! 4/24/12: replace the 1.563 with x**3, here and in nCAR routine.
              IF ( diag_id%qnidt_ihom +  diag_id%qni_ihom_col > 0 ) &
                    diag_4l(i,j,k,diag_pt%qnidt_ihom) =    &
                           dum*3._r8*dumc(i,k)/   &
                                (4._r8*3.14_r8*1.563e-14_r8*500._r8)/deltat
#endif
              end if
           end if

           ! remove any excess over-saturation, which is possible due to non-linearity when adding 
           ! together all microphysical processes
           !-----------------------------------------------------------------
           ! follow code similar to old CAM scheme

           qtmp=q(i,k)+qvlat(i,k)*deltat
           ttmp=t(i,k)+tlat(i,k)/cpp*deltat

           ! use rhw to allow ice supersaturation
#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (ttmp, p(i,k), qsn, q=q(i,k), &
                             esat = esn, es_over_liq = .true.)
#else
           call qsat_water(ttmp, p(i,k), esn, qsn)
           qsn = min(qsn,1._r8)
#endif

           if (qtmp > qsn .and. qsn > 0) then
              ! expression below is approximate since there may be ice deposition
              dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat
              ! add to output cme
              cmeout(i,k) = cmeout(i,k)+dum

              ! now add to tendencies, partition between liquid and ice based on temperature
#ifdef GFDL_COMPATIBLE_MICROP
        if (.not. clubb_active) then
          if (tiedtke_macrophysics) then
            if (ttmp > tmelt - 5._r8) then
                dum1 = 0.0_r8
                ssat_disposal(i,k) = 1._r8
              else if (ttmp < tmelt - 40._r8) then
                dum1 = 1.0_r8
                ssat_disposal(i,k) = 2._r8
              else
                dum1 = 0.0_r8
                ssat_disposal(i,k) = 1._r8
              end if
          else ! (tiedtke_macrophysics)

! for non-tiedtke, need to define how supersaturation removal affects
! particle number, if at all  ??
! for now, assume supersaturation removal does not lead to change in
! particle numbers
! use ice / liq partitioning as in original NCAR
              ssat_disposal(i,k) = 0._r8
              if (ttmp > tmelt-5._r8) then
                 dum1=0.0_r8
                 ! now add to tendencies, partition between liquid and ice based on te
                 !-------------------------------------------------------
              else if (ttmp < tmelt-35._r8) then
                 dum1=1.0_r8
              else
                 dum1=(tmelt-5._r8-ttmp)/30._r8
              end if
          endif ! (tiedtke_macrophysics)

      else  !clubb_active
#endif
           if (ttmp > 268.15_r8) then
              dum1=0.0_r8
! now add to tendencies, partition between liquid and ice based on te
           else if (ttmp < 238.15_r8) then
              dum1=1.0_r8
           else
              dum1=(268.15_r8-ttmp)/30._r8
           end if
#ifdef GFDL_COMPATIBLE_MICROP
       endif
#endif

              dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
                   *qsn/(cpp*rv*ttmp**2))/deltat
              qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
              ! for output
              qcrestot(i,k)=dum*(1._r8-dum1)
              qitend(i,k)=qitend(i,k)+dum*dum1
              qirestot(i,k)=dum*dum1
              qvlat(i,k)=qvlat(i,k)-dum
              ! for output
              qvres(i,k)=-dum
              tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
#ifdef GFDL_COMPATIBLE_MICROP
          else
              if (.not. clubb_active) &
            ssat_disposal(i,k) = 0._r8
#endif
          end if

        end if

        ! calculate effective radius for pass to radiation code
        !=========================================================
        ! if no cloud water, default value is 10 micron for droplets,
        ! 25 micron for cloud ice

        ! update cloud variables after instantaneous processes to get effective radius
        ! variables are in-cloud to calculate size dist parameters

        dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
        dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)

        ! hm add 6/2/11 switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k)=ncnst/rho(i,k)
        end if

        ! hm add 6/2/11 switch for specification of cloud ice number
        if (nicons) then
           dumni(i,k)=ninst/rho(i,k)
        end if

        ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1

        dumc(i,k)=min(dumc(i,k),5.e-3_r8)
        dumi(i,k)=min(dumi(i,k),5.e-3_r8)

        ! cloud ice effective radius
        !-----------------------------------------------------------------

        if (do_cldice) then
           if (dumi(i,k).ge.qsmall) then
#ifdef GFDL_COMPATIBLE_MICROP
           IF (diag_id%qnidt_size_adj +  diag_id%qni_size_adj_col > 0) &
                 diag_4l(i,j,k,diag_pt%qnidt_size_adj ) =  nitend(i,k)
#endif

              dum_2D(i,k) = dumni(i,k)
              call size_dist_param_ice(dumi(i,k), dumni(i,k), lami(i,k), n0i(i,k))

              if (dumni(i,k) /=dum_2D(i,k)) then
                 ! adjust number conc if needed to keep mean size in reasonable range
                 nitend(i,k)=(dumni(i,k)*icldm(i,k)-ni(i,k))/deltat
              end if

#ifdef GFDL_COMPATIBLE_MICROP
            IF (diag_id%qnidt_size_adj +  diag_id%qni_size_adj_col > 0) &
                 diag_4l(i,j,k,diag_pt%qnidt_size_adj ) =     &
                     nitend(i,k) - diag_4l(i,j,k,diag_pt%qnidt_size_adj )
#endif
              effi(i,k) = 1.5_r8/lami(i,k)*1.e6_r8

           else
              effi(i,k) = 25._r8
           end if

           ! ice effective diameter for david mitchell's optics
           deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8
        else
           ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
           ! radius has already been determined from the size distribution.
           effi(i,k) = re_ice(i,k) * 1.e6_r8      ! m -> um
           deffi(i,k)=effi(i,k) * 2._r8
        end if

        ! cloud droplet effective radius
        !-----------------------------------------------------------------
        if (dumc(i,k).ge.qsmall) then

#ifdef GFDL_COMPATIBLE_MICROP
           IF (diag_id%qndt_size_adj + diag_id%qn_size_adj_col  > 0) &
                 diag_4l(i,j,k,diag_pt%qndt_size_adj ) =  nctend(i,k)
#endif

           ! hm add 6/2/11 switch for specification of droplet and crystal number
           if (nccons) then
              ! make sure nc is consistence with the constant N by adjusting tendency, need
              ! to multiply by cloud fraction
              ! note that nctend may be further adjusted below if mean droplet size is
              ! out of bounds

              nctend(i,k)=(ncnst/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat

           end if           

           dum = dumnc(i,k)

           call size_dist_param_liq(dumc(i,k), dumnc(i,k), cdnl, rho(i,k), .true., &
                pgam(i,k), lamc(i,k))

           if (dum /= dumnc(i,k)) then
              ! adjust number conc if needed to keep mean size in reasonable range
              nctend(i,k)=(dumnc(i,k)*lcldm(i,k)-nc(i,k))/deltat
           end if
 
#ifdef GFDL_COMPATIBLE_MICROP
            IF (diag_id%qndt_size_adj + diag_id%qn_size_adj_col  > 0) &
                 diag_4l(i,j,k,diag_pt%qndt_size_adj ) =     &
                        nctend(i,k) - diag_4l(i,j,k,diag_pt%qndt_size_adj )
#endif

           effc(i,k) = gamma(pgam(i,k)+4._r8)/ &
                gamma(pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8
           !assign output fields for shape here
           lamcrad(i,k)=lamc(i,k)
           pgamrad(i,k)=pgam(i,k)


           ! recalculate effective radius for constant number, in order to separate
           ! first and second indirect effects
           !======================================
           ! assume constant number of 10^8 kg-1

           dumnc(i,k)=1.e8_r8

           ! Pass in "false" adjust flag to prevent number from being changed within
           ! size distribution subroutine.
           call size_dist_param_liq(dumc(i,k), dumnc(i,k), cdnl, rho(i,k), .false., &
                pgam(i,k), lamc(i,k))

           effc_fn(i,k) = gamma(pgam(i,k)+4._r8)/ &
                gamma(pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8

        else
           effc(i,k) = 10._r8
           lamcrad(i,k)=0._r8
           pgamrad(i,k)=0._r8
           effc_fn(i,k) = 10._r8
        end if

     end do ! vertical k loop

     do k=1,nlev
        ! if updated q (after microphysics) is zero, then ensure updated n is also zero
        !=================================================================================

#ifdef GFDL_COMPATIBLE_MICROP
          IF (diag_id%qndt_fill2  + diag_id%qn_fill2_col > 0) &
               diag_4l(i,j,k,diag_pt%qndt_fill2 ) =  nctend(i,k)
          IF (diag_id%qnidt_fill2 +  diag_id%qni_fill2_col > 0) &
               diag_4l(i,j,k,diag_pt%qnidt_fill2 ) =  nitend(i,k)
#endif
        if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat
        if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat

#ifdef GFDL_COMPATIBLE_MICROP
          IF (diag_id%qndt_fill2  + diag_id%qn_fill2_col > 0) &
               diag_4l(i,j,k,diag_pt%qndt_fill2 ) =     &
                          nctend(i,k) - diag_4l(i,j,k,diag_pt%qndt_fill2 )
          IF (diag_id%qnidt_fill2 +  diag_id%qni_fill2_col > 0) &
               diag_4l(i,j,k,diag_pt%qnidt_fill2 ) =     &
                         nitend(i,k) - diag_4l(i,j,k,diag_pt%qnidt_fill2 )
#endif
     end do

  end do sed_col_loop! i loop

  ! DO STUFF FOR OUTPUT:
  !==================================================

  ! qc and qi are only used for output calculations past here,
  ! so add qctend and qitend back in one more time
  qc = qc + qctend*deltat
  qi = qi + qitend*deltat

  ! averaging for snow and rain number and diameter
  !--------------------------------------------------

  ! drout2/dsout2:
  ! diameter of rain and snow
  ! dsout:
  ! scaled diameter of snow (passed to radiation in CAM)
  ! reff_rain/reff_snow:
  ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual

  where (qrout .gt. 1.e-7_r8 &
       .and. nrout.gt.0._r8)
     qrout2 = qrout * cldmax
     nrout2 = nrout * cldmax
     ! The avg_diameter call does the actual calculation; other diameter
     ! outputs are just drout2 times constants.
     drout2 = avg_diameter(qrout, nrout, rho, rhow) * cldmax
     freqr = cldmax

#ifdef GFDL_COMPATIBLE_MICROP
!RSH output instantaneous values for use in radiation rather than a 
!    cldfrac-weighted value which would require normalization by freqr 
!    (which also would need to be output).
     reff_rain=1.5_r8*drout2*1.e6_r8/ cldmax 
!    if( clubb_active ) then
!RSH routine now returns effective RADIUS, not DIAMETER. Conversion to
!   diameter now must be made in calling routine.
!       lsc_rain_size(i,k) = max(  60.0_r8, lsc_rain_size(i,k) )
!       lsc_rain_size(i,k) = min(1500.0_r8, lsc_rain_size(i,k) )
!       reff_rain          = max(  30.0_r8, reff_rain          )
!       reff_rain          = min( 750.0_r8, reff_rain          )
!    endif
#else
     reff_rain=1.5_r8*drout2*1.e6_r8 
#endif
  elsewhere
     qrout2 = 0._r8
     nrout2 = 0._r8
     drout2 = 0._r8
     freqr = 0._r8
#ifdef GFDL_COMPATIBLE_MICROP
     reff_rain = 50._r8
#else
     reff_rain = 0._r8
#endif
  end where
#ifdef GFDL_COMPATIBLE_MICROP
     if( clubb_active ) then
!RSH routine now returns effective RADIUS, not DIAMETER. Conversion to
!   diameter now must be made in calling routine.
!       lsc_rain_size(i,k) = max(  60.0_r8, lsc_rain_size(i,k) )
!       lsc_rain_size(i,k) = min(1500.0_r8, lsc_rain_size(i,k) )
        reff_rain          = max(  30.0_r8, reff_rain          )
        reff_rain          = min( 750.0_r8, reff_rain          )
     endif
#endif

  where (qsout .gt. 1.e-7_r8 &
       .and. nsout.gt.0._r8)
     qsout2 = qsout * cldmax
     nsout2 = nsout * cldmax
     ! The avg_diameter call does the actual calculation; other diameter
     ! outputs are just dsout2 times constants.
     dsout2 = avg_diameter(qsout, nsout, rho, rhosn)
     freqs = cldmax

     dsout=3._r8*rhosn/917._r8*dsout2

     dsout2 = dsout2 * cldmax

#ifdef GFDL_COMPATIBLE_MICROP
!RSH routine now returns effective RADIUS, not DIAMETER. Conversion to
!   diameter now must be made in calling routine.
!RSH output instantaneous values for use in radiation rather than a 
!    cldfrac-weighted value which would require normalization by freqs 
!    (which also would need to be output). Explore this issue later.
     reff_snow=1.5_r8*dsout2*1.e6_r8/cldmax
#else
     reff_snow=1.5_r8*dsout2*1.e6_r8
#endif
  elsewhere
     dsout  = 0._r8
     qsout2 = 0._r8
     nsout2 = 0._r8
     dsout2 = 0._r8
     freqs  = 0._r8
     reff_snow=0._r8
  end where

  ! analytic radar reflectivity
  !--------------------------------------------------
  ! formulas from Matthew Shupe, NOAA/CERES
  ! *****note: radar reflectivity is local (in-precip average)
  ! units of mm^6/m^3

  do i = 1,mgncol
     do k=1,nlev
        if (qc(i,k).ge.qsmall) then
           dum=(qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
                /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/cldmax(i,k)
        else
           dum=0._r8
        end if
        if (qi(i,k).ge.qsmall) then
           dum1=(qi(i,k)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/cldmax(i,k)
        else 
           dum1=0._r8
        end if

        if (qsout(i,k).ge.qsmall) then
           dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
        end if

        refl(i,k)=dum+dum1

        ! add rain rate, but for 37 GHz formulation instead of 94 GHz
        ! formula approximated from data of Matrasov (2007)
        ! rainrt is the rain rate in mm/hr
        ! reflectivity (dum) is in DBz

        if (rainrt(i,k).ge.0.001_r8) then
           dum=log10(rainrt(i,k)**6._r8)+16._r8

           ! convert from DBz to mm^6/m^3

           dum = 10._r8**(dum/10._r8)
        else
           ! don't include rain rate in R calculation for values less than 0.001 mm/hr
           dum=0._r8
        end if

        ! add to refl

        refl(i,k)=refl(i,k)+dum

        !output reflectivity in Z.
        areflz(i,k)=refl(i,k) * cldmax(i,k)

        ! convert back to DBz 

        if (refl(i,k).gt.minrefl) then 
           refl(i,k)=10._r8*log10(refl(i,k))
        else
           refl(i,k)=-9999._r8
        end if

        !set averaging flag
        if (refl(i,k).gt.mindbz) then 
           arefl(i,k)=refl(i,k) * cldmax(i,k)
           frefl(i,k)=cldmax(i,k)
        else
           arefl(i,k)=0._r8
           areflz(i,k)=0._r8
           frefl(i,k)=0._r8
        end if

        ! bound cloudsat reflectivity

        csrfl(i,k)=min(csmax,refl(i,k))

        !set averaging flag
        if (csrfl(i,k).gt.csmin) then 
           acsrfl(i,k)=refl(i,k) * cldmax(i,k)
           fcsrfl(i,k)=cldmax(i,k)
        else
           acsrfl(i,k)=0._r8
           fcsrfl(i,k)=0._r8
        end if

     end do
  end do

  !redefine fice here....
  dum_2D = qsout + qrout + qc + qi
  dumi = qsout + qi
  where (dumi .gt. qsmall .and. dum_2D .gt. qsmall)
     nfice=min(dumi/dum_2D,1._r8)
  elsewhere
     nfice=0._r8
  end where

  ! Unpack all outputs

  ! Avoid zero/near-zero division.
  qcsinksum_rate1ord = qcsinksum_rate1ord/max(qcsum_rate1ord,1.0e-30_r8)

  call unpack_array(qcsinksum_rate1ord, mgcols, top_lev, 0._r8, rate1ord_cw2pr_st)

  if(allocated(rndst))  deallocate(rndst)    !--->h1g, 2017-02-02
  if(allocated(nacon))  deallocate(nacon)    !--->h1g, 2017-02-02

#ifdef GFDL_COMPATIBLE_MICROP
  call unpack_array(tlat, mgcols, top_lev, 0._r8, tlato)
  call unpack_array(qvlat, mgcols, top_lev, 0._r8, qvlato)

  call unpack_array(qctend, mgcols, top_lev, 0._r8, qctendo)
  call unpack_array(qitend, mgcols, top_lev, 0._r8, qitendo)

  ! Note that where there is no water, we set nctend and nitend to remove number
  ! concentration as well.
  call unpack_array(nctend, mgcols, top_lev, 0._r8      , nctendo)
  if (do_cldice) then
     call unpack_array(nitend, mgcols, top_lev, 0._r8      , nitendo)
  else
     call unpack_array(nitend, mgcols, top_lev, 0._r8, nitendo)
  end if
  call unpack_array(prect, mgcols, 0._r8, precto)
  call unpack_array(preci, mgcols, 0._r8, precio)
  call unpack_array(qrout, mgcols, top_lev, 0._r8, qrouto)
  call unpack_array(qsout, mgcols, top_lev, 0._r8, qsouto)
  call unpack_array(reff_rain, mgcols, top_lev, 0._r8, reff_raino)
  call unpack_array(reff_snow, mgcols, top_lev, 0._r8, reff_snowo)
  call unpack_array(rflx, mgcols, top_lev, 0._r8, rflxo)
  call unpack_array(sflx, mgcols, top_lev, 0._r8, sflxo)
  call unpack_array(ssat_disposal, mgcols, top_lev, 0._r8, ssat_disposalout)
#endif
#ifndef GFDL_COMPATIBLE_MICROP
  call unpack_array(tlat, mgcols, top_lev, 0._r8, tlato)
  call unpack_array(qvlat, mgcols, top_lev, 0._r8, qvlato)

  call unpack_array(qctend, mgcols, top_lev, 0._r8, qctendo)
  call unpack_array(qitend, mgcols, top_lev, 0._r8, qitendo)

  ! Note that where there is no water, we set nctend and nitend to remove number
  ! concentration as well.
  call unpack_array(nctend, mgcols, top_lev, -ncn/deltat, nctendo)
  if (do_cldice) then
     call unpack_array(nitend, mgcols, top_lev, -nin/deltat, nitendo)
  else
     call unpack_array(nitend, mgcols, top_lev, 0._r8, nitendo)
  end if

  call unpack_array(effc, mgcols, top_lev, 10._r8, effco)
  call unpack_array(effc_fn, mgcols, top_lev, 10._r8, effco_fn)
  call unpack_array(effi, mgcols, top_lev, 25._r8, effio)

  call unpack_array(prect, mgcols, 0._r8, precto)
  call unpack_array(preci, mgcols, 0._r8, precio)

  call unpack_array(nevapr, mgcols, top_lev, 0._r8, nevapro)
  call unpack_array(evapsnow, mgcols, top_lev, 0._r8, evapsnowo)
  call unpack_array(prain, mgcols, top_lev, 0._r8, praino)
  call unpack_array(prodsnow, mgcols, top_lev, 0._r8, prodsnowo)
  call unpack_array(cmeout, mgcols, top_lev, 0._r8, cmeouto)

  call unpack_array(lamcrad, mgcols, top_lev, 0._r8, lamcrado)
  call unpack_array(pgamrad, mgcols, top_lev, 0._r8, pgamrado)
  call unpack_array(deffi, mgcols, top_lev, 0._r8, deffio)

  call unpack_array(qsout, mgcols, top_lev, 0._r8, qsouto)
  call unpack_array(qsout2, mgcols, top_lev, 0._r8, qsouto2)
  call unpack_array(nsout, mgcols, top_lev, 0._r8, nsouto)
  call unpack_array(nsout2, mgcols, top_lev, 0._r8, nsouto2)
  call unpack_array(dsout, mgcols, top_lev, 0._r8, dsouto)
  call unpack_array(dsout2, mgcols, top_lev, 0._r8, dsouto2)

  call unpack_array(qrout, mgcols, top_lev, 0._r8, qrouto)
  call unpack_array(qrout2, mgcols, top_lev, 0._r8, qrouto2)
  call unpack_array(nrout, mgcols, top_lev, 0._r8, nrouto)
  call unpack_array(nrout2, mgcols, top_lev, 0._r8, nrouto2)
  call unpack_array(drout2, mgcols, top_lev, 0._r8, drouto2)

  call unpack_array(reff_rain, mgcols, top_lev, 0._r8, reff_raino)
  call unpack_array(reff_snow, mgcols, top_lev, 0._r8, reff_snowo)

  call unpack_array(freqs, mgcols, top_lev, 0._r8, freqso)
  call unpack_array(freqr, mgcols, top_lev, 0._r8, freqro)

  call unpack_array(rflx, mgcols, top_lev, 0._r8, rflxo)
  call unpack_array(sflx, mgcols, top_lev, 0._r8, sflxo)

  call unpack_array(qcsevap, mgcols, top_lev, 0._r8, qcsevapo)
  call unpack_array(qisevap, mgcols, top_lev, 0._r8, qisevapo)
  call unpack_array(qvres, mgcols, top_lev, 0._r8, qvreso)
  call unpack_array(cmeitot, mgcols, top_lev, 0._r8, cmeiout)
  call unpack_array(vtrmc, mgcols, top_lev, 0._r8, vtrmco)
  call unpack_array(vtrmi, mgcols, top_lev, 0._r8, vtrmio)
  call unpack_array(qcsedten, mgcols, top_lev, 0._r8, qcsedteno)
  call unpack_array(qisedten, mgcols, top_lev, 0._r8, qisedteno)

  call unpack_array(pratot, mgcols,top_lev, 0._r8, prao)
  call unpack_array(prctot, mgcols,top_lev, 0._r8, prco)
  call unpack_array(mnuccctot, mgcols,top_lev, 0._r8, mnuccco)
  call unpack_array(mnuccttot, mgcols,top_lev, 0._r8, mnuccto)
  call unpack_array(msacwitot, mgcols,top_lev, 0._r8, msacwio)
  call unpack_array(psacwstot, mgcols,top_lev, 0._r8, psacwso)
  call unpack_array(bergstot, mgcols,top_lev, 0._r8, bergso)
  call unpack_array(bergtot, mgcols,top_lev, 0._r8, bergo)
  call unpack_array(melttot, mgcols,top_lev, 0._r8, melto)
  call unpack_array(homotot, mgcols,top_lev, 0._r8, homoo)
  call unpack_array(qcrestot, mgcols,top_lev, 0._r8, qcreso)
  call unpack_array(prcitot, mgcols,top_lev, 0._r8, prcio)
  call unpack_array(praitot, mgcols,top_lev, 0._r8, praio)
  call unpack_array(qirestot, mgcols,top_lev, 0._r8, qireso)
  call unpack_array(mnuccrtot, mgcols,top_lev, 0._r8, mnuccro)
  call unpack_array(pracstot, mgcols,top_lev, 0._r8, pracso)
  call unpack_array(mnuccdtot, mgcols,top_lev, 0._r8, mnuccdo)
  call unpack_array(meltsdttot, mgcols,top_lev, 0._r8, meltsdto)
  call unpack_array(frzrdttot, mgcols,top_lev, 0._r8, frzrdto)

  call unpack_array(refl, mgcols, top_lev, -9999._r8, reflo)
  call unpack_array(arefl, mgcols, top_lev, 0._r8, areflo)
  call unpack_array(areflz, mgcols, top_lev, 0._r8, areflzo)
  call unpack_array(frefl, mgcols, top_lev, 0._r8, freflo)
  call unpack_array(csrfl, mgcols, top_lev, -9999._r8, csrflo)
  call unpack_array(acsrfl, mgcols, top_lev, 0._r8, acsrflo)
  call unpack_array(fcsrfl, mgcols, top_lev, 0._r8, fcsrflo)

  call unpack_array(rercld, mgcols, top_lev, 0._r8, rercldo)

  call unpack_array(nfice, mgcols, top_lev, 0._r8, nficeo)

  call unpack_array(ncai, mgcols, top_lev, 0._r8, ncaio)
  call unpack_array(ncal, mgcols, top_lev, 0._r8, ncalo)
#endif

#ifdef GFDL_COMPATIBLE_MICROP
! diagnostics for water tendencies
! water  vapor specific humicity
      if  (diag_id%rain_evap + diag_id%rain_evap_col > 0)  &
              diag_4l(:,j,:,diag_pt%rain_evap)  = -preo( : , : )
      if  (diag_id%qdt_rain_evap > 0)  &
              diag_4l(:,j,:,diag_pt%qdt_rain_evap)  = -preo( : , : )
      if (diag_id%qdt_cond   > 0) &
              diag_4l(:,j,:,diag_pt%qdt_cond)  = -cmelo(:,:)
      if  (diag_id%qdt_snow_sublim > 0 .or.   &
                                     diag_id%q_snow_sublim_col > 0 )  &
              diag_4l(:,j,:,diag_pt%qdt_snow_sublim )  =  - prdso( : , : )
      if  (diag_id%qdt_deposition > 0)  &
              diag_4l(:,j,:,diag_pt%qdt_deposition )  = -cmeitot( : , : )
      if  (diag_id%qdt_sedi_ice2vapor> 0)  &
              diag_4l(:,j,:,diag_pt%qdt_sedi_ice2vapor) = qisevap( : , : )
      if  (diag_id%qdt_sedi_liquid2vapor> 0)  &
            diag_4l(:,j,:,diag_pt%qdt_sedi_liquid2vapor) = qcsevap( : , : )
      if  (diag_id%qdt_super_sat_rm > 0)  &
              diag_4l(:,j,:,diag_pt%qdt_super_sat_rm) = qvres( : , : )

! cloud liquid water
      if (diag_id%qldt_accr  + diag_id%ql_accr_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_accr)  = - pratot(:,:)
      if (diag_id%qldt_auto  + diag_id%ql_auto_col > 0)&
              diag_4l(:,j,:,diag_pt%qldt_auto)  = -prctot(:,:)
      if (diag_id%qldt_freez2 + diag_id%ql_freez2_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_freez2) =   &
                                          -(mnuccctot(:,:) + mnuccttot(:,:) )
      sum_freeze2(:,:) =  mnuccctot(:,:) + mnuccttot(:,:)
      if (diag_id%qldt_accrs  + diag_id%ql_accrs_col > 0) & 
              diag_4l(:,j,:,diag_pt%qldt_accrs)  = -psacwstot(:,:) 
      sum_rime(:,:) =  psacwstot(:,:)
      if (diag_id%qldt_HM_splinter + diag_id%ql_HM_splinter_col > 0)&
              diag_4l(:,j,:,diag_pt%qldt_HM_splinter)  = -msacwitot(:,:)
      sum_splinter(:,:) =  msacwitot(:,:)
      if (diag_id%qldt_bergs + diag_id%ql_bergs_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_bergs)  = -bergstot(:,:)
      sum_bergs(:,:) =  bergstot (:,:)
      if (diag_id%qidt_dep + diag_id%qi_dep_col > 0)    &
              diag_4l(:,j,:,diag_pt%qidt_dep)  = max(cmeitot(:,:),0._r8) 
      sum_cond(:,:) = max(cmeitot(:,:),0._r8) 
      if (diag_id%qidt_subl + diag_id%qi_subl_col > 0)  &
              diag_4l(:,j,:,diag_pt%qidt_subl)  =     &
                                         -max(-1._r8*cmeitot(:,:),0._r8) 
      if (diag_id%qldt_cond + diag_id%ql_cond_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_cond)  =  max(cmelo(:,:), 0._r8)
      if (diag_id%qldt_evap  + diag_id%ql_evap_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_evap)  =   &
                                           - max(-1._r8*cmelo(:,:),0._r8)
      if (diag_id%qldt_eros + diag_id%ql_eros_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_eros)  =   eroslo   (:,:)
      if (diag_id%qdt_eros_l                       > 0) &
              diag_4l(:,j,:,diag_pt%qdt_eros_l)  =  -eroslo   (:,:)
      if (diag_id%qidt_eros + diag_id%qi_eros_col > 0) &
              diag_4l(:,j,:,diag_pt%qidt_eros)  =  erosio   (:,:)
      if (diag_id%qdt_eros_i                       > 0) &
              diag_4l(:,j,:,diag_pt%qdt_eros_i)  = -erosio   (:,:)
      if (diag_id%qldt_berg + diag_id%ql_berg_col > 0) &
              diag_4l(:,j,:,diag_pt%qldt_berg)  =  -bergtot(:,:)
      sum_berg(:,:) =  bergtot(:,:)
      IF ( diag_id%qldt_sedi  + diag_id%ql_sedi_col > 0 ) &
              diag_4l(:,j,1:nlev,diag_pt%qldt_sedi) = qcsedten(:,1:nlev)
      IF ( diag_id%liq_adj  + diag_id%liq_adj_col > 0 ) &
              diag_4l(:,j,:,diag_pt%liq_adj) = qcrestot(:,:)

! cloud ice water
      if (diag_id%qidt_auto + diag_id%qi_auto_col > 0) &
             diag_4l(:,j,:,diag_pt%qidt_auto) = -prcitot(:,:)
      if (diag_id%qidt_accr  + diag_id%qi_accr_col > 0) &
             diag_4l(:,j,:,diag_pt%qidt_accr) = -praitot(:,:)
      IF ( diag_id%qidt_fall  + diag_id%qi_fall_col > 0 ) &
              diag_4l(:,j,1:nlev,diag_pt%qidt_fall) = qisedten(:,1:nlev)
      IF ( diag_id%ice_adj  +  diag_id%ice_adj_col > 0 ) &
              diag_4l(:,j,:,diag_pt%ice_adj) = qirestot(:,:)
      sum_ice_adj(:,:) = qirestot(:,:)

! ---> rain water mixing ratio
      if (diag_id%srfrain_accrs + diag_id%srfrain_accrs_col > 0)    &
             diag_4l(:,j,:,diag_pt%srfrain_accrs)  = -pracstot(:,:)
      if (diag_id%srfrain_freez + diag_id%srfrain_freez_col > 0)    &
             diag_4l(:,j,:,diag_pt%srfrain_freez)  = -mnuccrtot(:,:)

! ---> snow mixing ratio

! --->liquid droplet number
      if (diag_id%qndt_cond + diag_id%qn_cond_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_cond)  = npccno(:,:)  /real(iter)
      if (diag_id%qndt_freez + diag_id%qn_freez_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_freez)  = nnuccco(:,:)  /real(iter)
      if (diag_id%qndt_contact_frz > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_contact_frz)  =    &
                                                 nnuccto(:,:) /real(iter)
      if (diag_id%qndt_sacws + diag_id%qn_sacws_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_sacws)  = npsacwso(:,:) /real(iter)
      if (diag_id%qndt_evap + diag_id%qn_evap_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_evap)  = nsubco(:,:)  /real(iter)
      if (diag_id%qndt_eros + diag_id%qn_eros_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_eros)  = nerosco(:,:)  /real(iter)
      if (diag_id%qndt_pra + diag_id%qn_pra_col > 0)    &
             diag_4l(:,j,:,diag_pt%qndt_pra)  = nprao(:,:)  /real(iter)
      if (diag_id%qndt_auto + diag_id%qn_auto_col > 0)    &
              diag_4l(:,j,:,diag_pt%qndt_auto)  = nprc1o(:,:)  /real(iter)
      if ( diag_id%qndt_nucclim  + diag_id%qn_nucclim_col  > 0 ) &
            diag_4l(:,j,:,diag_pt%qndt_nucclim)  =     &
                                                nucclimo(:,:)  /real(iter)

! ---> ice number 
      if (diag_id%qnidt_nnuccd +  diag_id%qni_nnuccd_col > 0)    &
             diag_4l(:,j,:,diag_pt%qnidt_nnuccd)  =    &
                                                 nnuccdo(:,:)  /real(iter)
      if (diag_id%qnidt_nsacwi> 0)    &
             diag_4l(:,j,:,diag_pt%qnidt_nsacwi)  =     &
                                                 nsacwio(:,:)  /real(iter)
      if (diag_id%qnidt_nsubi  + diag_id%qni_nsubi_col  > 0)    &
             diag_4l(:,j,:,diag_pt%qnidt_nsubi)  = nsubio(:,:)  /real(iter)
      if (diag_id%qnidt_nerosi  + diag_id%qni_nerosi_col  > 0)    &
           diag_4l(:,j,:,diag_pt%qnidt_nerosi)  = nerosio(:,:)  /real(iter)
      if (diag_id%qnidt_nprci  + diag_id%qni_nprci_col  > 0)    &
             diag_4l(:,j,:,diag_pt%qnidt_nprci)  = nprcio(:,:)  /real(iter)
      if (diag_id%qnidt_nprai  + diag_id%qni_nprai_col  > 0)    &
             diag_4l(:,j,:,diag_pt%qnidt_nprai)  = npraio(:,:)  /real(iter)
      if (diag_id%qnidt_nucclim1 +  diag_id%qni_nucclim1_col > 0 ) &
             diag_4l(:,j,:,diag_pt%qnidt_nucclim1) =      &
                                              nucclim1io(:,:)  /real(iter)
!10/23/13: NOTE: STILL NEEDS CONVERSION !!!!
!RSH:
!   calculate fraction of total ice / snow creation that requires 
!   ice-forming nuclei
      do k=1,nlev
        do i=1,mgncol
          qldt_sum = sum_cond(i,k) + sum_rime(i,k) + sum_berg(i,k) + &
                     sum_ice_adj(i,k) + MAX(sum_bergs(i,k), 0.0) + &
                     sum_freeze(i,k) + sum_freeze2(i,k) + sum_splinter(i,k)
          if (ABS(qldt_sum) > 0.0            ) then
! ---> h1g, 2014-07-18, add option of including contact freeze in bergeron
           if( include_contact_freeze_in_berg ) then
            f_snow_berg_l(i,k) = (sum_berg(i,k) + sum_cond(i,k) +   &
                                  sum_ice_adj(i,k) +    &
                                  MAX( sum_bergs(i,k), 0.0) +     &
                                  sum_freeze (i,k) + sum_freeze2(i,k) )/qldt_sum
           else
            f_snow_berg_l(i,k) = (sum_berg(i,k) + sum_cond(i,k) +   &
                                sum_ice_adj(i,k) +    &
                                MAX( sum_bergs(i,k), 0.0) +     &
                                sum_freeze (i,k) + mnuccctot(i,k))/qldt_sum        
            endif
! <--- h1g, 2014-07-18
          else
            f_snow_berg_l(i,k) = 0._r8
          endif
        end do
      end do
      call unpack_array (f_snow_berg_l,mgcols, top_lev, 0._r8, f_snow_berg)
      call unpack_diag_array (diag_4l(:,j,:,0:), mgcols, j, top_lev, diag_4d(:,j,:,0:))

#endif
end subroutine micro_mg_tend

!========================================================================
!FORMULAS
!========================================================================

! Calculate correction due to latent heat for evaporation/sublimation
elemental function calc_ab(t, qv, xxl) result(ab)
  real(r8), intent(in) :: t     ! Temperature
  real(r8), intent(in) :: qv    ! Saturation vapor pressure
  real(r8), intent(in) :: xxl   ! Latent heat

  real(r8) :: ab

  real(r8) :: dqsdt

  dqsdt = xxl*qv / (rv * t**2)
  ab = 1._r8 + dqsdt*xxl/cpp

end function calc_ab

! get cloud droplet size distribution parameters
#ifdef GFDL_COMPATIBLE_MICROP
           subroutine size_dist_param_liq(qcic, ncic, cdnl, rho, nadjflag, pgam, lamc)
#else
 elemental subroutine size_dist_param_liq(qcic, ncic, cdnl, rho, nadjflag, pgam, lamc)
#endif

  real(r8), intent(in) :: qcic
  real(r8), intent(inout) :: ncic
  real(r8), intent(in) :: cdnl
  real(r8), intent(in) :: rho
  logical,  intent(in) :: nadjflag ! Whether to adjust number concentration to fall
                                   ! within certain bounds

  real(r8), intent(out) :: pgam
  real(r8), intent(out) :: lamc

  real(r8) :: dumgam1
  real(r8) :: dumgam2
  real(r8) :: lammin
  real(r8) :: lammax

  if (qcic > qsmall) then

     if (nadjflag) then
        ! add upper limit to in-cloud number concentration to prevent numerical error
        ncic=min(ncic,qcic*1.e20_r8)
        ! add lower limit to in-cloud number concentration
        ncic=max(ncic,cdnl/rho) ! sghan minimum in #/cm
     end if

     ! get pgam from fit to observations of martin et al. 1994

     pgam=0.0005714_r8*(ncic/1.e6_r8*rho)+0.2714_r8
     pgam=1._r8/(pgam**2)-1._r8
     pgam=max(pgam,2._r8)
     pgam=min(pgam,15._r8)

     ! calculate lamc
     dumgam1 = gamma(pgam+1._r8)
     dumgam2 = gamma(pgam+4._r8)

     lamc = (pi/6._r8*rhow*ncic*dumgam2/ &
          (qcic*dumgam1))**(1._r8/3._r8)

     ! lammin, 50 micron diameter max mean size

     lammin = (pgam+1._r8)/50.e-6_r8
     lammax = (pgam+1._r8)/2.e-6_r8

     if (lamc < lammin) then
        lamc = lammin

        if (nadjflag) then
           ncic = 6._r8 * lamc**3 * qcic * dumgam1/ &
                (pi * rhow * dumgam2)
        end if
     else if (lamc > lammax) then
        lamc = lammax

        if (nadjflag) then
           ncic = 6._r8 * lamc**3 * qcic * dumgam1/ &
                (pi * rhow * dumgam2)
        end if
     end if

  else
     ! pgam not calculated in this case, so set it to a value likely to cause an error
     ! if it is accidentally used
     ! (gamma function undefined for negative integers)
     pgam = -100._r8
     lamc = 0._r8
  end if

end subroutine size_dist_param_liq

! get ice size distribution parameters
elemental subroutine size_dist_param_ice(qiic, niic, lami, n0i)
  real(r8), intent(in) :: qiic
  real(r8), intent(inout) :: niic

  real(r8), intent(out) :: lami
  real(r8), intent(out) :: n0i


  ! local parameters
#ifdef GFDL_COMPATIBLE_MICROP
  real(r8)            :: lammini 
  real(r8)            :: lammaxi 
  real(r8)            :: ci
#else
  real(r8), parameter :: lammaxi = 1._r8/10.e-6_r8
  real(r8), parameter :: lammini = 1._r8/(2._r8*dcs)
  ! particle mass-diameter relationship
  ! currently we assume spherical particles for cloud ice/snow
  ! m = cD^d
  ! cloud ice mass-diameter relationship
  real(r8), parameter :: ci = rhoi*pi/6._r8
#endif

#ifdef GFDL_COMPATIBLE_MICROP
  lammini = 1._r8/(2._r8*dcs)
  lammaxi = 1._r8/min_diam_ice
  ci = rhoi*pi/6._r8
#endif

  if (qiic > qsmall) then

     ! add upper limit to in-cloud number concentration to prevent numerical error
     niic = min(niic, qiic * 1.e20_r8)

     lami = (cons1*ci*niic/qiic)**(1._r8/dsph)
     n0i = niic * lami
  
     ! check for slope
     ! adjust vars
     if (lami < lammini) then
        lami = lammini
        n0i = lami**(dsph+1._r8) * qiic/(ci*cons1)
        niic = n0i/lami
     else if (lami > lammaxi) then
        lami = lammaxi
        n0i = lami**(dsph+1._r8) * qiic/(ci*cons1)
        niic = n0i/lami
     end if
  else
     lami = 0._r8
     n0i  = 0._r8
  end if

end subroutine size_dist_param_ice

! get rain size distribution parameters
elemental subroutine size_dist_param_rain(qric, nric, lamr, n0r)
  real(r8), intent(in) :: qric
  real(r8), intent(inout) :: nric

  real(r8), intent(out) :: lamr
  real(r8), intent(out) :: n0r

  ! particle mass-diameter relationship
  ! currently we assume spherical particles for cloud ice/snow
  ! m = cD^d
  ! rain mass-diameter relationship
  ! Rain is hard-coded for spherical drops
  ! Therefore cr is rhow*pi, rather than
  ! using rhow*pi/6 and later multiplying
  ! by gamma(1+dsph) == 6.
  real(r8), parameter :: cr = rhow*pi

  ! local parameters
  real(r8), parameter :: lammaxr = 1._r8/20.e-6_r8
  real(r8), parameter :: lamminr = 1._r8/500.e-6_r8

  if (qric > qsmall) then

     lamr = (cr*nric/qric)**(1._r8/3._r8)
     n0r = nric * lamr
  
     ! check for slope
     ! adjust vars

     if (lamr < lamminr) then
        lamr = lamminr
        n0r = lamr**4 * qric/cr
        nric = n0r/lamr
     else if (lamr > lammaxr) then
        lamr = lammaxr
        n0r = lamr**4 * qric/cr
        nric = n0r/lamr
     end if

  else
     lamr = 0._r8
     n0r  = 0._r8
  end if

end subroutine size_dist_param_rain

! get snow size distribution parameters
elemental subroutine size_dist_param_snow(qsic, nsic, lams, n0s)
  real(r8), intent(in) :: qsic
  real(r8), intent(inout) :: nsic

  real(r8), intent(out) :: lams
  real(r8), intent(out) :: n0s

  ! particle mass-diameter relationship
  ! currently we assume spherical particles for cloud ice/snow
  ! m = cD^d
  ! cloud ice mass-diameter relationship
#ifdef GFDL_COMPATIBLE_MICROP
  real(r8)            :: cs 
#else
  real(r8), parameter :: cs = rhosn*pi/6._r8
#endif


#ifdef GFDL_COMPATIBLE_MICROP
  real(r8) :: lammaxs 
  real(r8) :: lammins
#else
  ! local parameters
  real(r8), parameter :: lammaxs = 1._r8/10.e-6_r8
  real(r8), parameter :: lammins = 1._r8/2000.e-6_r8
#endif


#ifdef GFDL_COMPATIBLE_MICROP
   cs = rhosn*pi/6._r8
   lammaxs = 1._r8/min_diam_ice
   lammins = 1._r8/2000.e-6_r8
#endif

  if (qsic > qsmall) then

     lams = (cons1*cs*nsic/qsic)**(1._r8/dsph)
     n0s = nsic * lams
  
     ! check for slope
     ! adjust vars
     if (lams < lammins) then
        lams = lammins
        n0s = lams**(dsph+1._r8) * qsic/(cs*cons1)
        nsic = n0s/lams
     else if (lams > lammaxs) then
        lams = lammaxs
        n0s = lams**(dsph+1._r8) * qsic/(cs*cons1)
        nsic = n0s/lams
     end if
  
  else
     lams = 0._r8
     n0s  = 0._r8
  end if

end subroutine size_dist_param_snow

real(r8) elemental function avg_diameter(q, n, rho_air, rho_sub)
  ! Finds the average diameter of particles given their density, and
  ! mass/number concentrations in the air.
  real(r8), intent(in) :: q         ! mass mixing ratio
  real(r8), intent(in) :: n         ! number concentration
  real(r8), intent(in) :: rho_air   ! local density of the air
  real(r8), intent(in) :: rho_sub   ! density of the particle substance

  avg_diameter = (pi * rho_sub * n/(q*rho_air))**(-1._r8/3._r8)

end function avg_diameter

#ifdef GFDL_COMPATIBLE_MICROP
real(r8)           function var_coef(relvar, a)
#else
real(r8) elemental function var_coef(relvar, a)
#endif
  ! Finds a coefficient for process rates based on the relative variance
  ! of cloud water.
  real(r8), intent(in) :: relvar
  real(r8), intent(in) :: a

  var_coef = gamma(relvar + a) / (gamma(relvar) * relvar**a)

end function var_coef

!========================================================================
!MICROPHYSICAL PROCESS CALCULATIONS
!========================================================================

!========================================================================
! Initial ice deposition and sublimation loop.
! Run before the main loop

#ifdef GFDL_COMPATIBLE_MICROP
elemental subroutine ice_deposition_sublimation_init(clubb_active, deltat, t, q, qc, qi, ni, &
                                                lcldm, icldm, naai, rho, dv, &
                                                esl, esi, qvl, qvi, relhum, &
                                                berg, cmei, cmel, &
                                                tiedtke_macrophysics,  &
                                                dqcdt, dqidt)
#else
elemental subroutine ice_deposition_sublimation_init(deltat, t, q, qc, qi, ni, &
                                                lcldm, icldm, naai, rho, dv, &
                                                esl, esi, qvl, qvi, relhum, &
                                                berg, cmei)
#endif

  ! Inputs
  logical, intent(in)  :: clubb_active
  real(r8), intent(in) :: deltat

  real(r8), intent(in) :: t
  real(r8), intent(in) :: q

  real(r8), intent(in) :: qc
  real(r8), intent(in) :: qi
  real(r8), intent(in) :: ni

  real(r8), intent(in) :: lcldm
  real(r8), intent(in) :: icldm

  real(r8), intent(in) :: naai
  real(r8), intent(in) :: rho
  real(r8), intent(in) :: dv

  real(r8), intent(in) :: esl
  real(r8), intent(in) :: esi
  real(r8), intent(in) :: qvl
  real(r8), intent(in) :: qvi
  real(r8), intent(in) :: relhum

#ifdef GFDL_COMPATIBLE_MICROP
  real(r8), intent(in) :: dqcdt, dqidt
  logical, intent(in) :: tiedtke_macrophysics
  real(r8), intent(out) :: cmel
  logical  :: do_berg1
  logical  :: limit_berg 
  real(r8) :: berg_lim 
  real(r8) :: dum2
#endif
  ! Outputs

  real(r8), intent(out) :: berg
  real(r8), intent(out) :: cmei

  ! Internal variables

  real(r8) :: qiic
  real(r8) :: niic

  real(r8) :: prd               ! provisional deposition rate of cloud ice at water sat 
  real(r8) :: bergtsf           ! bergeron timescale to remove all liquid
  real(r8) :: rhin              ! modified RH for vapor deposition

  real(r8) :: epsi              ! 1/sat relaxation timescale for ice

  real(r8) :: ab
  real(r8) :: lami
  real(r8) :: n0i
  real(r8) :: dum

  ! initialize bergeron process to zero
  berg = 0._r8

  ! Initialize CME components
  cmei = 0._r8
  
  if (t < icenuct) then
     ! provisional nucleation rate
     dum = max((naai - ni/icldm)/deltat*icldm,0._r8)
  else
     dum = 0._r8
  end if

  ! get in-cloud qi and ni after nucleation
#ifdef GFDL_COMPATIBLE_MICROP
  if (tiedtke_macrophysics) then
    qiic = qi/icldm 
  else
    qiic = (qi + dum*deltat*mi0)/icldm
  endif
#else
  qiic = (qi + dum*deltat*mi0)/icldm
#endif

  niic = (ni + dum*deltat)/icldm

  ! hm add 6/2/11 switch for specification of cloud ice number
  if (nicons) niic = ninst/rho

#ifdef GFDL_COMPATIBLE_MICROP
!  define the large-scale cloud water and ice condensation/evaporation from
!  values passed into routine, and their sum.
          limit_berg = .false.
          berg_lim = 1.0e-6_r8
          cmel = dqcdt
          cmei = dqidt
          dum2 = cmel + cmei

!  if bergeron process to be active for any non-zero condensate, set flag
!  so indicating.
!  current setting is limit_berg = .false. (controlled by nml)
          IF ( .NOT.  limit_berg ) THEN
            if (dum2 .ge. 0._r8) then
              do_berg1 = .true.
            else
              do_berg1 = .false.
           end if
          ELSE
! GFDL has option to not allow bergeron process when cloud ice is less than
! berg_lim, even if have positive condensate. set flag appropriately. 
!           if (dum2 .ge. 0._r8 .and. qinew .gt. berg_lim ) then
            if (dum2 .ge. 0._r8 .and. qiic*icldm .gt. berg_lim ) then
              do_berg1 = .true.
            else
              do_berg1 = .false.
            end if
          END If

          if (do_berg1 .or. clubb_active) THEN
#endif

 !ICE DEPOSITION:
  !=============================================
  !if ice exists
  if (t < tmelt .and. qi >= qsmall) then

     ab = calc_ab(t, qvi, xxls)

     ! get ice size distribution parameters
     call size_dist_param_ice(qiic, niic, lami, n0i)

     epsi = 2._r8*pi*n0i*rho*Dv/(lami*lami)
              
     !if liquid exists
#ifdef GFDL_COMPATIBLE_MICROP
!RSHCLUBB:
! is CLUBB entering dqcdt = 0. always?? ; if not there needs to be 
! an clubb_active loop here
     if (qc + dqcdt*deltat >= qsmall) then
#else
     if (qc >= qsmall) then
#endif
                 
        ! calculate Bergeron process
        berg = epsi*(qvl-qvi)/ab * Bergeron_fac  ! h1g, 2014-07-20
                 
        ! multiply by cloud fraction
        berg = berg*min(icldm,lcldm)
                 
        ! Must be positive
        if (berg <= 0._r8) then
           berg = 0._r8
        else
           !BERGERON LIMITING WHEN ALL LIQUID DEPLETED IN 1 TIMESTEP
           !-------------------------------------------------------------
                 
#ifdef GFDL_COMPATIBLE_MICROP
      if( .not. clubb_active ) then
           bergtsf = (qc/deltat+dqcdt)/berg  ! bergeron time scale (fraction of timestep)
      else
#endif
           bergtsf = (qc/berg) / deltat ! bergeron time scale (fraction of timestep)
#ifdef GFDL_COMPATIBLE_MICROP
      endif
#endif
           if (bergtsf < 1._r8) then

#ifdef GFDL_COMPATIBLE_MICROP
      if( .not. clubb_active ) then
              berg = dqcdt + qc/deltat
      else
#endif
              berg = qc/deltat
#ifdef GFDL_COMPATIBLE_MICROP
      endif
#endif
                    
            
#ifdef GFDL_COMPATIBLE_MICROP
!  As per Marc, the  following is inconsistent with the Tiedtke asumption 
!  of in-cloud RH = 1., so it is excluded in the Tiedtke case.
              if (.not. tiedtke_macrophysics) then
#endif
                 
              rhin  = (1.0_r8 + relhum) / 2._r8
              !assume RH for frac of step w/ no liq is 1/2 way btwn cldy & cell-ave RH.

              if ((rhin*esl/esi) > 1._r8) then !if ice saturated (but all liquid evap'd)
                 prd = epsi*(rhin*qvl-qvi)/ab

                 ! multiply by cloud fraction assuming liquid/ice maximum overlap
                 prd = prd*min(icldm,lcldm)
   
                 ! add to cmei
                 cmei = cmei + (prd * (1._r8- bergtsf))

              end if ! rhin 
                    
#ifdef GFDL_COMPATIBLE_MICROP
            endif ! (not tiedtke_macrophysics)
#endif
         end if  ! bergtsf
                    
        end if   ! berg <=0
        end if   ! qc >= qsmall

#ifdef GFDL_COMPATIBLE_MICROP
      if (.not. tiedtke_macrophysics) then
#endif
     !Ice deposition in frac of cell with no liquid
     !-------------------------------------------------------------
     ! store liquid cloud fraction in 'dum'
     if (qc >= qsmall) then
        dum = lcldm
     else
        ! for case of no liquid, need to set liquid cloud fraction to zero
        dum = 0._r8
     end if

     if (icldm > dum) then

        ! set RH to grid-mean value for pure ice cloud
        rhin = relhum
                 
        if ((rhin*esl/esi) > 1._r8) then !if rh over ice>ice saturation.
                    
           prd = epsi*(rhin*qvl-qvi)/ab
                    
           ! multiply by relevant cloud fraction for pure ice cloud
           ! assuming maximum overlap of liquid/ice
           prd = prd*(icldm-dum) !apply to ice-only part of cld.
           cmei = cmei + prd
                    
        end if ! rhin
     end if ! qc or icldm > lcldm

     !if grid-mean is ice saturated & qi formed in non-liq cld part, 
     !limit ice formation to avoid mean becoming undersaturated.
     !-------------------------------------------------------------
     if(cmei > 0.0_r8 .and. (relhum*esl/esi) > 1._r8 ) &
          ! max berg is val which removes all ice supersaturation from vapor phase. 
          cmei=min(cmei,(q-qvl*esi/esl)/ab/deltat)
#ifdef GFDL_COMPATIBLE_MICROP
    endif  ! .not. tiedtke_macrophysics
#endif
#ifdef GFDL_COMPATIBLE_MICROP
! Marc includes a restriction on berg at T < -40C
                if (t < tmelt  - 40._r8 .and. ( .not. clubb_active) ) then
                  berg = 0._r8
                endif
#endif

  end if ! end ice exists and t < tmelt

#ifdef GFDL_COMPATIBLE_MICROP
  endif  ! (do_berg1 .or. clubb_active)
#endif
  !ICE SUBLIMATION:
#ifdef GFDL_COMPATIBLE_MICROP
  if (.not. tiedtke_macrophysics) then
#endif
  !=========================================
  !If ice-subsaturated and ice exists:

  if ((relhum*esl/esi) < 1._r8 .and. qiic >= qsmall ) then

     ab = calc_ab(t, qvi, xxls)

     ! get ice size distribution parameters
     call size_dist_param_ice(qiic, niic, lami, n0i)

     epsi = 2._r8*pi*n0i*rho*Dv/(lami*lami)

     ! modify for ice fraction below
     prd = epsi*(relhum*qvl-qvi)/ab * icldm
     cmei=min(prd,0._r8)
           
  endif !subsaturated and ice exists

  ! sublimation should not exceed available ice
  cmei = max(cmei, -qi/deltat)
        
  ! sublimation should not increase grid mean rhi above 1.0 
  if(cmei < 0.0_r8 .and. (relhum*esl/esi) < 1._r8 ) &
       cmei=min(0._r8,max(cmei,(q-qvl*esi/esl)/ab/deltat))
#ifdef GFDL_COMPATIBLE_MICROP
  endif  ! (not. tiedtke_macrophysics) 
#endif

  ! limit cmei due for roundoff error
  cmei = cmei*omsm
#ifdef GFDL_COMPATIBLE_MICROP
  if (.not. clubb_active) &
      cmel = cmel*omsm
#endif

end subroutine ice_deposition_sublimation_init

!========================================================================
! autoconversion of cloud liquid water to rain
! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
! minimum qc of 1 x 10^-8 prevents floating point error

#ifdef GFDL_COMPATIBLE_MICROP
          subroutine kk2000_liq_autoconversion(qcic, ncic, rho, relvar, &
#else
elemental subroutine kk2000_liq_autoconversion(qcic, ncic, rho, relvar, &
#endif
                        prc, nprc, nprc1)
  
  real(r8), intent(in) :: qcic
  real(r8), intent(in) :: ncic
  real(r8), intent(in) :: rho
  real(r8), intent(in) :: relvar

  real(r8), intent(out) :: prc
  real(r8), intent(out) :: nprc
  real(r8), intent(out) :: nprc1

  real(r8) :: prc_coef
  real(r8) :: nprc_denom

  ! subcolumn modifications change coefficients
  if (microp_uniform) then
     prc_coef = 1._r8
     nprc_denom = (4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
  else
     prc_coef = var_coef(relvar, 2.47_r8)
     nprc_denom = cons22
  end if

  q_cri = 4._r8/3._r8*pi*rhow*(r_cri**3) * ncic
  q_cri = max (icsmall, q_cri)
  
  if (qcic .ge. q_cri) then

     ! nprc is increase in rain number conc due to autoconversion
     ! nprc1 is decrease in cloud droplet conc due to autoconversion
     
     ! assume exponential sub-grid distribution of qc, resulting in additional
     ! factor related to qcvar below
     ! hm switch for sub-columns, don't include sub-grid qc
     
      if ( use_K2013_auto ) then
         prc = prc_coef * &
            7.98e10_r8 * qcic**4.22_r8 * (ncic/1.e6_r8*rho)**(-3.01_r8)
      elseif ( use_RK_auto ) then
        prc =  &
          32681. * qcic**2.33_r8 * (ncic*rho*rhow)**(-0.33_r8)
      else
     prc = prc_coef * &
          1350._r8 * qcic**2.47_r8 * (ncic/1.e6_r8*rho)**(-1.79_r8)
      endif
     nprc = prc/nprc_denom
     nprc1 = prc/(qcic/ncic)

  else
     prc=0._r8
     nprc=0._r8
     nprc1=0._r8
  end if

end subroutine kk2000_liq_autoconversion

!========================================================================
! Autoconversion of cloud ice to snow
! similar to Ferrier (1994)

elemental subroutine ice_autoconversion(t, qiic, lami, n0i, prci, nprci)

  real(r8), intent(in) :: t
  real(r8), intent(in) :: qiic
  real(r8), intent(in) :: lami
  real(r8), intent(in) :: n0i
  
  real(r8), intent(out) :: prci
  real(r8), intent(out) :: nprci

  if (t .le. tmelt .and.qiic.ge.qsmall) then

     ! note: assumes autoconversion timescale of 180 sec

     nprci = n0i/(lami*180._r8)*exp(-lami*dcs)

     prci = pi*rhoi*n0i/(6._r8*180._r8)* &
          (cons23/lami+3._r8*cons24/lami**2+ &
          6._r8*dcs/lami**3+6._r8/lami**4)*exp(-lami*dcs)          

  else
     prci=0._r8
     nprci=0._r8
  end if

end subroutine ice_autoconversion

! immersion freezing (Bigg, 1953)
!===================================

#ifdef GFDL_COMPATIBLE_MICROP
          subroutine immersion_freezing(t, pgam, lamc, cdist1, qcic, &
#else
elemental subroutine immersion_freezing(t, pgam, lamc, cdist1, qcic, &
#endif
                        relvar, mnuccc, nnuccc)

  ! Temperature
  real(r8), intent(in) :: t

  ! Cloud droplet size distribution parameters
  real(r8), intent(in) :: pgam
  real(r8), intent(in) :: lamc
  real(r8), intent(in) :: cdist1

  ! MMR of in-cloud liquid water
  real(r8), intent(in) :: qcic
  
  ! Relative variance of cloud water
  real(r8), intent(in) :: relvar

  ! Output tendencies
  real(r8), intent(out) :: mnuccc ! MMR
  real(r8), intent(out) :: nnuccc ! Number
  
  ! Coefficients that will be omitted for sub-columns
  real(r8) :: dum, dum1


  if (microp_uniform) then
     dum = 1._r8
     dum1 = 1._r8
  else
     dum = var_coef(relvar, 2._r8)
     dum1 = var_coef(relvar, 1._r8)
  end if

#ifdef GFDL_COMPATIBLE_MICROP
  if (qcic >= qsmall .and. t < tmelt-4._r8) then
#else
  if (qcic >= qsmall .and. t < 269.15_r8) then
#endif

     mnuccc = dum * &
          pi*pi/36._r8*rhow* &
          cdist1*gamma(7._r8+pgam)* &
          bimm*(exp(aimm*(tmelt - t))-1._r8)/lamc**3/lamc**3

     nnuccc = dum1 * &
          pi/6._r8*cdist1*gamma(pgam+4._r8) &
          *bimm*(exp(aimm*(tmelt - t))-1._r8)/lamc**3

  else
     mnuccc = 0._r8
     nnuccc = 0._r8
  end if ! qcic > qsmall and t < 4 deg C

end subroutine immersion_freezing

! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
!===================================================================
! dust size and number in multiple bins are read in from companion routine

!pure subroutine contact_freezing (t, p, rndst, nacon, pgam, lamc, cdist1, qcic, &
     subroutine contact_freezing (t, p, rndst, nacon, pgam, lamc, cdist1, qcic, &
     relvar, mnucct, nnucct)
  
  real(r8), intent(in) :: t(:)            ! Temperature
  real(r8), intent(in) :: p(:)            ! Pressure
  real(r8), intent(in) :: rndst(:,:)      ! Radius (for multiple dust bins)
  real(r8), intent(in) :: nacon(:,:)      ! Number (for multiple dust bins)

  ! Size distribution parameters for cloud droplets
  real(r8), intent(in) :: pgam(:)
  real(r8), intent(in) :: lamc(:)
  real(r8), intent(in) :: cdist1(:)

  ! MMR of in-cloud liquid water
  real(r8), intent(in) :: qcic(:)
  
  ! Relative cloud water variance
  real(r8), intent(in) :: relvar(:)

  ! Output tendencies
  real(r8), intent(out) :: mnucct(:) ! MMR
  real(r8), intent(out) :: nnucct(:) ! Number

  real(r8) :: tcnt                  ! scaled relative temperature
  real(r8) :: viscosity             ! temperature-specific viscosity (kg/m/s)
  real(r8) :: mfp                   ! temperature-specific mean free path (m)

  ! Dimension these according to number of dust bins, inferred from rndst size
  real(r8) :: nslip(size(rndst,2))  ! slip correction factors
  real(r8) :: ndfaer(size(rndst,2)) ! aerosol diffusivities (m^2/sec)

  ! Coefficients not used for subcolumns
  real(r8) :: dum, dum1

  integer  :: i

  ! subcolumns

  do i = 1,size(t)

     if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then

        if (microp_uniform) then
           dum = 1._r8
           dum1 = 1._r8
        else
           dum = var_coef(relvar(i), 4._r8/3._r8)
           dum1 = var_coef(relvar(i), 1._r8/3._r8)
        endif

        tcnt=(270.16_r8-t(i))**1.3_r8
        viscosity = 1.8e-5_r8*(t(i)/298.0_r8)**0.85_r8    ! Viscosity (kg/m/s)
        mfp = 2.0_r8*viscosity/ &                         ! Mean free path (m)
                     (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) ))

        ! Note that these two are vectors.
        nslip = 1.0_r8+(mfp/rndst(i,:))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,:)/mfp))))! Slip correction factor
        ndfaer = 1.381e-23_r8*t(i)*nslip/(6._r8*pi*viscosity*rndst(i,:))  ! aerosol diffusivity (m2/s)

        mnucct(i) = dum *  &
             dot_product(ndfaer,nacon(i,:)*tcnt)*pi*pi/3._r8*rhow* &
             cdist1(i)*gamma(pgam(i)+5._r8)/lamc(i)**4

        nnucct(i) =  dum1 *  &
             dot_product(ndfaer,nacon(i,:)*tcnt)*2._r8*pi*  &
             cdist1(i)*gamma(pgam(i)+2._r8)/lamc(i)

     else

        mnucct(i)=0._r8
        nnucct(i)=0._r8

     end if ! qcic > qsmall and t < 4 deg C
  end do

end subroutine contact_freezing

! snow self-aggregation from passarelli, 1978, used by reisner, 1998
!===================================================================
! this is hard-wired for bs = 0.4 for now
! ignore self-collection of cloud ice

elemental subroutine snow_self_aggregation(t, rho, asn, qsic, nsic, nsagg)

  real(r8), intent(in) :: t    ! Temperature
  real(r8), intent(in) :: rho  ! Density
  real(r8), intent(in) :: asn  ! fall speed parameter for snow

  ! In-cloud snow
  real(r8), intent(in) :: qsic ! MMR
  real(r8), intent(in) :: nsic ! Number

  ! Output number tendency
  real(r8), intent(out) :: nsagg

  if (qsic >= qsmall .and. t <= tmelt) then
     nsagg = -1108._r8*asn*eii* &
          pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)* &
          rho**((2._r8+bs)/3._r8)*qsic**((2._r8+bs)/3._r8)* &
          (nsic*rho)**((4._r8-bs)/3._r8) /(4._r8*720._r8*rho)
  else
     nsagg=0._r8
  end if

end subroutine snow_self_aggregation

! accretion of cloud droplets onto snow/graupel
!===================================================================
! here use continuous collection equation with
! simple gravitational collection kernel
! ignore collisions between droplets/cloud ice
! since minimum size ice particle for accretion is 50 - 150 micron

elemental subroutine accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, &
     pgam, lamc, lams, n0s, psacws, npsacws)

  real(r8), intent(in) :: t   ! Temperature
  real(r8), intent(in) :: rho ! Density
  real(r8), intent(in) :: asn ! Fallspeed parameter (snow)
  real(r8), intent(in) :: uns ! Current fallspeed   (snow)
  real(r8), intent(in) :: mu  ! Viscosity

  ! In-cloud liquid water
  real(r8), intent(in) :: qcic ! MMR
  real(r8), intent(in) :: ncic ! Number

  ! In-cloud snow
  real(r8), intent(in) :: qsic ! MMR

  ! Cloud droplet size parameters
  real(r8), intent(in) :: pgam
  real(r8), intent(in) :: lamc

  ! Snow size parameters
  real(r8), intent(in) :: lams
  real(r8), intent(in) :: n0s

  ! Output tendencies
  real(r8), intent(out) :: psacws  ! Mass mixing ratio
  real(r8), intent(out) :: npsacws ! Number concentration

  real(r8) :: dc0 ! Provisional mean droplet size
  real(r8) :: dum
  real(r8) :: eci ! collection efficiency for riming of snow by droplets

  ! ignore collision of snow with droplets above freezing

  if (qsic >= qsmall .and. t <= tmelt .and. qcic >= qsmall) then

     ! put in size dependent collection efficiency
     ! mean diameter of snow is area-weighted, since
     ! accretion is function of crystal geometric area
     ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)

     dc0 = (pgam+1._r8)/lamc
     dum = dc0*dc0*uns*rhow/(9._r8*mu*(1._r8/lams))
     eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))

     eci = max(eci,0._r8)
     eci = min(eci,1._r8)

     ! no impact of sub-grid distribution of qc since psacws
     ! is linear in qc

     psacws = pi/4._r8*asn*qcic*rho*n0s*eci*cons11 / lams**(bs+3._r8)
     npsacws = pi/4._r8*asn*ncic*rho*n0s*eci*cons11 / lams**(bs+3._r8)
  else
     psacws = 0._r8
     npsacws = 0._r8
  end if

end subroutine accrete_cloud_water_snow

! add secondary ice production due to accretion of droplets by snow 
!===================================================================
! (Hallet-Mossop process) (from Cotton et al., 1986)

elemental subroutine secondary_ice_production(t, psacws, msacwi, nsacwi)
  real(r8), intent(in) :: t ! Temperature

  ! Accretion of cloud water to snow tendencies
  real(r8), intent(inout) :: psacws ! MMR

  ! Output (ice) tendencies
  real(r8), intent(out) :: msacwi ! MMR
  real(r8), intent(out) :: nsacwi ! Number

  if((t < 270.16_r8) .and. (t >= 268.16_r8)) then
     nsacwi = 3.5e8_r8*(270.16_r8-t)/2.0_r8*psacws
     msacwi = min(nsacwi*mi0, psacws)
  else if((t < 268.16_r8) .and. (t >= 265.16_r8)) then
     nsacwi = 3.5e8_r8*(t-265.16_r8)/3.0_r8*psacws
     msacwi = min(nsacwi*mi0, psacws)
  else
     nsacwi = 0.0_r8
     msacwi = 0.0_r8
  endif

  psacws = max(0.0_r8,psacws - nsacwi*mi0)

end subroutine secondary_ice_production

! accretion of rain water by snow
!===================================================================
! formula from ikawa and saito, 1991, used by reisner et al., 1998

elemental subroutine accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, &
     lamr, n0r, lams, n0s, pracs, npracs )

  real(r8), intent(in) :: t   ! Temperature
  real(r8), intent(in) :: rho ! Density

  ! Fallspeeds
  ! mass-weighted
  real(r8), intent(in) :: umr ! rain
  real(r8), intent(in) :: ums ! snow
  ! number-weighted
  real(r8), intent(in) :: unr ! rain
  real(r8), intent(in) :: uns ! snow

  ! In cloud MMRs
  real(r8), intent(in) :: qric ! rain
  real(r8), intent(in) :: qsic ! snow

  ! Size distribution parameters
  ! rain
  real(r8), intent(in) :: lamr
  real(r8), intent(in) :: n0r
  ! snow
  real(r8), intent(in) :: lams
  real(r8), intent(in) :: n0s

  ! Output tendencies
  real(r8), intent(out) :: pracs  ! MMR
  real(r8), intent(out) :: npracs ! Number

  ! Collection efficiency for accretion of rain by snow
  real(r8), parameter :: ecr = 1.0_r8

  if (qric >= icsmall .and. qsic >= icsmall .and. t <= tmelt) then

     pracs = pi*pi*ecr*(((1.2_r8*umr-0.95_r8*ums)**2 + &
          0.08_r8*ums*umr)**0.5_r8 *  &
          rhow * rho * n0r * n0s * &
          (5._r8/(lamr**6 * lams)+ &
          2._r8/(lamr**5 * lams**2)+ &
          0.5_r8/(lamr**4 * lams**3)))

     npracs = pi/2._r8*rho*ecr* (1.7_r8*(unr-uns)**2 + &
          0.3_r8*unr*uns)**0.5_r8 * &
          n0r*n0s* &
          (1._r8/(lamr**3 * lams)+ &
          1._r8/(lamr**2 * lams**2)+ &
          1._r8/(lamr * lams**3))

  else
     pracs = 0._r8
     npracs = 0._r8
  end if

end subroutine accrete_rain_snow

! heterogeneous freezing of rain drops
!===================================================================
! follows from Bigg (1953)

elemental subroutine heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr)

  real(r8), intent(in) :: t    ! Temperature

  ! In-cloud rain
  real(r8), intent(in) :: qric ! MMR
  real(r8), intent(in) :: nric ! Number
  real(r8), intent(in) :: lamr ! size parameter

  ! Output tendencies
  real(r8), intent(out) :: mnuccr ! MMR
  real(r8), intent(out) :: nnuccr ! Number

#ifdef GFDL_COMPATIBLE_MICROP
  if (t < tmelt-4._r8 .and. qric >= qsmall) then
#else
  if (t < 269.15_r8 .and. qric >= qsmall) then
#endif

     ! Division by lamr**3 twice is old workaround to avoid overflow.
     ! Probably no longer necessary
     mnuccr = 20._r8*pi*pi*rhow*nric*bimm* &
          (exp(aimm*(tmelt - t))-1._r8)/lamr**3 &
          /lamr**3

     nnuccr = pi*nric*bimm* &
          (exp(aimm*(tmelt - t))-1._r8)/lamr**3
  else
     mnuccr = 0._r8
     nnuccr = 0._r8
  end if
end subroutine heterogeneous_rain_freezing

! accretion of cloud liquid water by rain
!===================================================================
! formula from Khrouditnov and Kogan (2000)
! gravitational collection kernel, droplet fall speed neglected

#ifdef GFDL_COMPATIBLE_MICROP
          subroutine accrete_cloud_water_rain(qric, qcic, ncic, rho_air, &
                                              umr, &
#else
elemental subroutine accrete_cloud_water_rain(qric, qcic, ncic, rho_air, &
                                              umr, &
#endif
                     relvar, accre_enhan, pra, npra)

  ! In-cloud rain
  real(r8), intent(in) :: qric ! MMR

  ! Cloud droplets
  real(r8), intent(in) :: qcic ! MMR
  real(r8), intent(in) :: ncic ! Number
  
  ! SGS variability
  real(r8), intent(in) :: relvar
  real(r8), intent(in) :: accre_enhan

  ! air density
  real(r8), intent(in) :: rho_air
  ! Fallspeeds
  ! mass-weighted
  real(r8), intent(in) :: umr ! rain

  ! Output tendencies
  real(r8), intent(out) :: pra  ! MMR
  real(r8), intent(out) :: npra ! Number

  ! Coefficient that varies for subcolumns
  real(r8) :: pra_coef

  real(r8) :: rad_liq

  rad_liq = (0.75*qcic/ncic /pi /rhow)**(1.0/3.0)
  rad_liq = rad_liq*1.e6


  if (microp_uniform) then
     pra_coef = 1._r8
  else
     pra_coef = accre_enhan * var_coef(relvar, 1.15_r8)
  end if

  if (qric >= qsmall .and. qcic >= qsmall) then

     
     ! include sub-grid distribution of cloud water
    if ( use_K2013_accr ) then
     ! Kogan, JAS 2013, PP1423-1436
     pra = pra_coef * 8.51_r8*(qcic**1.05_r8)*(qric**0.98_r8)
    elseif ( use_RK_accr ) then
     pra =  65.772565* (rad_liq**2) / (rad_liq**2 + 20.5)* &
           ((qric*rho_air*umr/rhow)**(7.0/9.0)) * qcic
    else
     pra = pra_coef * 67._r8*(qcic*qric)**1.15_r8
    endif

     npra = pra/(qcic/ncic)

  else
     pra = 0._r8
     npra = 0._r8
  end if
end subroutine accrete_cloud_water_rain

! Self-collection of rain drops
!===================================================================
! from Beheng(1994)

elemental subroutine self_collection_rain(rho, qric, nric, nragg)

  real(r8), intent(in) :: rho  ! Air density

  ! Rain
  real(r8), intent(in) :: qric ! MMR
  real(r8), intent(in) :: nric ! Number

  ! Output number tendency
  real(r8), intent(out) :: nragg

  if (qric >= qsmall) then
     nragg = -8._r8*nric*qric*rho
  else
     nragg = 0._r8
  end if

end subroutine self_collection_rain

! Accretion of cloud ice by snow
!===================================================================
! For this calculation, it is assumed that the Vs >> Vi
! and Ds >> Di for continuous collection

elemental subroutine accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, &
     lams, n0s, prai, nprai)

  real(r8), intent(in) :: t    ! Temperature
  real(r8), intent(in) :: rho  ! Density

  real(r8), intent(in) :: asn  ! Snow fallspeed parameter

  ! Cloud ice
  real(r8), intent(in) :: qiic ! MMR
  real(r8), intent(in) :: niic ! Number

  real(r8), intent(in) :: qsic ! Snow MMR

  ! Snow size parameters
  real(r8), intent(in) :: lams
  real(r8), intent(in) :: n0s  

  ! Output tendencies
  real(r8), intent(out) :: prai  ! MMR
  real(r8), intent(out) :: nprai ! Number

  if (qsic >= qsmall .and. qiic >= qsmall .and. t <= tmelt) then

     prai = pi/4._r8 * asn * qiic * rho * n0s * eii * cons11/ &
          lams**(bs+3._r8)

     nprai = pi/4._r8 * asn * niic * rho * n0s * eii * cons11/ &
          lams**(bs+3._r8)
  else
     prai = 0._r8
     nprai = 0._r8
  end if

end subroutine accrete_cloud_ice_snow

! calculate evaporation/sublimation of rain and snow
!===================================================================
! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
! in-cloud condensation/deposition of rain and snow is neglected
! except for transfer of cloud water to snow through bergeron process

#ifdef GFDL_COMPATIBLE_MICROP
 subroutine evaporate_sublimate_precip(deltat, t, p, rho, dv, mu, sc, q, qvl, qvi, &
     lcldm, cldmax, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, cmei, D_eros_l, D_eros_i, cmel, pre, prds, clubb_active)
#else
 elemental subroutine evaporate_sublimate_precip(deltat, t, p, rho, dv, mu, sc, q, qvl, qvi, &
     lcldm, cldmax, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, cmei, pre, prds)
#endif

  real(r8), intent(in) :: deltat ! timestep

  real(r8), intent(in) :: t    ! temperature
  real(r8), intent(in) :: p    ! pressure
  real(r8), intent(in) :: rho  ! air density
  real(r8), intent(in) :: dv   ! water vapor diffusivity
  real(r8), intent(in) :: mu   ! viscosity
  real(r8), intent(in) :: sc   ! schmidt number
  real(r8), intent(in) :: q    ! humidity
  real(r8), intent(in) :: qvl  ! saturation humidity (water)
  real(r8), intent(in) :: qvi  ! saturation humidity (ice)
  real(r8), intent(in) :: lcldm  ! liquid cloud fraction
  real(r8), intent(in) :: cldmax ! precipitation fraction (maximum overlap)

  ! fallspeed parameters
  real(r8), intent(in) :: arn  ! rain
  real(r8), intent(in) :: asn  ! snow

  ! In-cloud MMRs
  real(r8), intent(in) :: qcic ! cloud liquid
  real(r8), intent(in) :: qiic ! cloud ice
  real(r8), intent(in) :: qric ! rain
  real(r8), intent(in) :: qsic ! snow

  ! Size parameters
  ! rain
  real(r8), intent(in) :: lamr
  real(r8), intent(in) :: n0r
  ! snow
  real(r8), intent(in) :: lams
  real(r8), intent(in) :: n0s

  ! cloud ice sublimation/deposition tendency
  real(r8), intent(in) :: cmei
#ifdef GFDL_COMPATIBLE_MICROP
  real(r8), intent(in) :: D_eros_l
  real(r8), intent(in) :: D_eros_i
  real(r8), intent(in) :: cmel
  logical,  intent(in) :: clubb_active
#endif

  ! Output tendencies
  real(r8), intent(out) :: pre
  real(r8), intent(out) :: prds

  ! checking for RH after rain evap
  real(r8) :: esn    ! saturation pressure
  real(r8) :: qsn    ! saturation humidity

  real(r8) :: qclr   ! water vapor mixing ratio in clear air
  real(r8) :: ab     ! correction to account for latent heat
  real(r8) :: eps    ! 1/ sat relaxation timescale

  ! Temps/dummies
  real(r8) :: qtmp
  real(r8) :: ttmp

  real(r8) :: dum, dum1

  ! set temporary cloud fraction to zero if cloud water + ice is very small
  ! this will ensure that evaporation/sublimation of precip occurs over
  ! entire grid cell, since min cloud fraction is specified otherwise
#ifdef GFDL_COMPATIBLE_MICROP
  if (qcic+qiic < qsmall) then
#else
  if (qcic+qiic < 1.e-6_r8) then
#endif
     dum = 0._r8
  else
     dum = lcldm
  end if

  ! only calculate if there is some precip fraction > cloud fraction

  if (cldmax > dum) then

     ! calculate q for out-of-cloud region
     qsn = min(qvl,1._r8)
     qclr=(q-dum*qsn)/(1._r8-dum)

     ! evaporation of rain
     if (qric.ge.qsmall) then

        ab = calc_ab(t, qvl, xxlv)
        eps = 2._r8*pi*n0r*rho*Dv* &
             (f1r/(lamr*lamr)+ &
             f2r*(arn*rho/mu)**0.5_r8* &
             sc**(1._r8/3._r8)*cons13/ &
             (lamr**(5._r8/2._r8+br/2._r8)))

        pre = eps*(qclr-qvl)/ab

        ! only evaporate in out-of-cloud region
        ! and distribute across cldmax
        pre=min(pre*(cldmax-dum),0._r8)
        pre=pre/cldmax
     else
        pre = 0._r8
     end if

     ! sublimation of snow
     if (qsic.ge.qsmall) then
        ab = calc_ab(t, qvi, xxls)
        eps = 2._r8*pi*n0s*rho*Dv* &
             (f1s/(lams*lams)+ &
             f2s*(asn*rho/mu)**0.5_r8* &
             sc**(1._r8/3._r8)*cons14/ &
             (lams**(5._r8/2._r8+bs/2._r8)))
        prds = eps*(qclr-qvi)/ab

        ! only sublimate in out-of-cloud region and distribute over cldmax
        prds=min(prds*(cldmax-dum),0._r8)
        prds=prds/cldmax
     else
        prds = 0._r8
     end if

     ! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
     ! get updated RH at end of time step based on cloud water/ice condensation/evap

#ifdef GFDL_COMPATIBLE_MICROP
       if (.not. clubb_active) then
              qtmp = q -     &
                       (D_eros_l + D_eros_i + cmel +   &
                        cmei + (pre + prds)*cldmax)*deltat
              ttmp = t +(  &
                       (D_eros_l + cmel + pre*cldmax)* xxlv + &
                       (D_eros_i + cmei + prds*cldmax)* xxls)*  &
                                                           deltat/cpp
       else
#endif
               qtmp=q-(cmei+(pre+prds)*cldmax)*deltat
               ttmp=t+((pre*cldmax)*xxlv+ &
                               (cmei+prds*cldmax)*xxls)*deltat/cpp
#ifdef GFDL_COMPATIBLE_MICROP
       endif
#endif

     !limit range of temperatures!
#ifdef GFDL_COMPATIBLE_MICROP
     ttmp = max(lowest_temp_for_sublimation, min(ttmp, 323._r8))
#else
     ttmp=max(180._r8,min(ttmp,323._r8))
#endif

     ! use rhw to allow ice supersaturation
#ifdef GFDL_COMPATIBLE_MICROP
     call compute_qs (ttmp, p, qsn, q = q, &
                       esat = esn, es_over_liq = .true.)
#else
     call qsat_water(ttmp, p, esn, qsn)
     qsn=min(qsn,1._r8)
#endif

     ! modify precip evaporation rate if q > qsat
     if (qtmp.gt.qsn) then
        if (pre+prds.lt.-1.e-20_r8) then
           dum1=pre/(pre+prds)
           ! recalculate q and t after cloud water cond but without precip evap
#ifdef GFDL_COMPATIBLE_MICROP
           if (.not. clubb_active) then
                  qtmp = q -    &
                               (D_eros_l + D_eros_i +   &
                                         cmel + cmei)*deltat
                  ttmp = t+    &
                            ((D_eros_l + cmel)*xxlv + &
                               (D_eros_i + cmei)*xxls)*deltat/cpp
           else
#endif
             qtmp=q-(cmei)*deltat
             ttmp=t+(cmei*xxls)*deltat/cpp
#ifdef GFDL_COMPATIBLE_MICROP
           endif
#endif
           ! use rhw to allow ice supersaturation
#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (ttmp, p, qsn, q = q, &
                       esat = esn, es_over_liq = .true.)
#else
           call qsat_water(ttmp, p, esn, qsn)
           qsn=min(qsn,1._r8)
#endif

           dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2))
           dum=min(dum,0._r8)

           ! modify rates if needed, divide by cldmax to get local (in-precip) value
           pre=dum*dum1/deltat/cldmax

           ! do separately using RHI for prds....
           ! use rhi to allow ice supersaturation
#ifdef GFDL_COMPATIBLE_MICROP
           call compute_qs (ttmp, p, qsn, q = q,   &
                             esat = esn, es_over_liq_and_ice = .true.)
#else
           call qsat_ice(ttmp, p, esn, qsn)
           qsn=min(qsn,1._r8)
#endif

           dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2))
           dum=min(dum,0._r8)

           ! modify rates if needed, divide by cldmax to get local (in-precip) value
           prds=dum*(1._r8-dum1)/deltat/cldmax
        end if
     end if

  else
     prds = 0._r8
     pre = 0._r8
  end if

end subroutine evaporate_sublimate_precip

! bergeron process - evaporation of droplets and deposition onto snow
!===================================================================

elemental subroutine bergeron_process(t, rho, dv, mu, sc, qvl, qvi, asn, &
     qcic, qsic, lams, n0s, bergs)

  real(r8), intent(in) :: t    ! temperature
  real(r8), intent(in) :: rho  ! air density
  real(r8), intent(in) :: dv   ! water vapor diffusivity
  real(r8), intent(in) :: mu   ! viscosity
  real(r8), intent(in) :: sc   ! schmidt number
  real(r8), intent(in) :: qvl  ! saturation humidity (water)
  real(r8), intent(in) :: qvi  ! saturation humidity (ice)

  ! fallspeed parameter for snow
  real(r8), intent(in) :: asn

  ! In-cloud MMRs
  real(r8), intent(in) :: qcic ! cloud liquid
  real(r8), intent(in) :: qsic ! snow

  ! Size parameters for snow
  real(r8), intent(in) :: lams
  real(r8), intent(in) :: n0s

  ! Output tendencies
  real(r8), intent(out) :: bergs

  real(r8) :: ab     ! correction to account for latent heat
  real(r8) :: eps    ! 1/ sat relaxation timescale

  if (qsic >= qsmall.and. qcic >= qsmall .and. t < tmelt) then
     ab = calc_ab(t, qvi, xxls)
     eps = 2._r8*pi*n0s*rho*Dv* &
          (f1s/(lams*lams)+ &
          f2s*(asn*rho/mu)**0.5_r8* &
          sc**(1._r8/3._r8)*cons14/ &
          (lams**(5._r8/2._r8+bs/2._r8)))
     bergs = eps*(qvl-qvi)/ab * Bergeron_fac  ! h1g, 2014-07-20
  else
     bergs = 0._r8
  end if

end subroutine bergeron_process

!========================================================================
!OUTPUT CALCULATIONS
!========================================================================

#ifdef GFDL_COMPATIBLE_MICROP
          subroutine calc_rercld(lamr, n0r, lamc, cdist1, pgam, dumr, qcic, &
#else
elemental subroutine calc_rercld(lamr, n0r, lamc, cdist1, pgam, dumr, qcic, &
#endif
     arcld, rercld)
  real(r8), intent(in) :: lamr          ! rain size parameter (slope)
  real(r8), intent(in) :: n0r           ! rain size parameter (intercept)
  real(r8), intent(in) :: lamc          ! size distribution parameter (slope)
  real(r8), intent(in) :: cdist1        ! for droplet freezing
  real(r8), intent(in) :: pgam          ! droplet size parameter
  real(r8), intent(in) :: dumr          ! in-cloud rain mass mixing ratio
  real(r8), intent(in) :: qcic          ! in-cloud cloud liquid

  integer,  intent(inout) :: arcld      ! number of substeps rercld has been through
  real(r8), intent(inout) :: rercld     ! effective radius calculation for rain + cloud

  ! combined size of precip & cloud drops
  real(r8) :: Atmp

  ! Rain drops
  if (lamr > 0._r8) then
     Atmp = n0r * pi / (2._r8 * lamr**3._r8)
  else
     Atmp = 0._r8
  end if

  ! Add cloud drops
  if (lamc > 0._r8) then
     Atmp = Atmp + cdist1 * pi * gamma(pgam+3._r8)/(4._r8 * lamc**2._r8)
  end if

  if (Atmp > 0._r8) then
     rercld = rercld + 3._r8 *(dumr + qcic) / (4._r8 * rhow * Atmp)
     arcld = arcld+1
  end if

end subroutine calc_rercld

!========================================================================
!UTILITIES
!========================================================================

#ifdef GFDL_COMPATIBLE_MICROP
pure subroutine micro_mg_get_cols(ncol, nlev, top_lev, qvn, qcn, qin, &
                   mgncol, mgcols, clubb_active)
#else
pure subroutine micro_mg_get_cols(ncol, nlev, top_lev, qcn, qin, &
     mgncol, mgcols)
#endif

  ! Determines which columns microphysics should operate over by
  ! checking for non-zero cloud water/ice.

  integer, intent(in) :: ncol      ! Number of columns with meaningful data
  integer, intent(in) :: nlev      ! Number of levels to use
  integer, intent(in) :: top_lev   ! Top level for microphysics

  real(r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg)
  real(r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg)
#ifdef GFDL_COMPATIBLE_MICROP
  real(r8), intent(in) :: qvn(:,:) ! cloud water mixing ratio (kg/kg)
  logical,  intent(in) :: clubb_active
#endif

  integer, intent(out) :: mgncol   ! Number of columns MG will use
  integer, allocatable, intent(out) :: mgcols(:) ! column indices

  integer :: lev_offset  ! top_lev - 1 (defined here for consistency)
  logical :: ltrue(ncol) ! store tests for each column

  integer :: i, ii ! column indices

  if (allocated(mgcols)) deallocate(mgcols)

  lev_offset = top_lev - 1

  ! Using "any" along dimension 2 collapses across levels, but
  ! not columns, so we know if water is present at any level
  ! in each column.

  ltrue = any(qcn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
  ltrue = ltrue .or. any(qin(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
#ifdef GFDL_COMPATIBLE_MICROP
 if (.not. clubb_active) then
  ltrue = ltrue .and.  .not. any(qvn(:ncol,top_lev:(nlev+lev_offset)) < -1.e-9_r8, 2)
  ltrue = ltrue .and.  .not. any(qvn(:ncol,top_lev:(nlev+lev_offset)) +    &
                                 qcn(:ncol,top_lev:(nlev+lev_offset)) +    &
                                 qin(:ncol,top_lev:(nlev+lev_offset) ) &
                                               < -1.e-9_r8, 2)
 endif ! (.not. clubb_active)
#endif

  ! Scan for true values to get a usable list of indices.
  
  mgncol = count(ltrue)
  allocate(mgcols(mgncol))
  i = 0
  do ii = 1,ncol
     if (ltrue(ii)) then
        i = i + 1
        mgcols(i) = ii
     end if
  end do

end subroutine micro_mg_get_cols

pure function interp_to_mid(orig_val, weights) result(new_val)
  ! Linear interpolation, here used to move from interfaces to midlevel
  real(r8), intent(in) :: orig_val(:,:)
  real(r8), intent(in) :: weights(:,:)

  ! New value is slightly smaller than the old value
  real(r8) :: new_val(size(orig_val,1),size(orig_val,2)-1)

  new_val = orig_val(:,:size(new_val,2))
  new_val = new_val + weights * (orig_val(:,2:) - new_val)

end function interp_to_mid

! Subroutines to pack arrays into smaller, contiguous pieces
!========================================================================
! Rank 1 array of reals, columns only
pure subroutine pack_array_1Dr8(old_array, cols, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:)   ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include

  ! Output
  real(r8), intent(out) :: new_array(:)
  
  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = old_array(cols)
  end if

end subroutine pack_array_1Dr8

! Rank 2 array of reals, columns and levels
pure subroutine pack_array_2Dr8(old_array, cols, top_lev, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  integer,  intent(in)  :: top_lev        ! First level to use

  ! Output
  real(r8), intent(out) :: new_array(:,:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = old_array(cols, top_lev:)
  end if

end subroutine pack_array_2Dr8

! Rank 3 array of reals, assume last index is extra
pure subroutine pack_array_3Dr8(old_array, cols, top_lev, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)          ! List of columns to include
  integer,  intent(in)  :: top_lev          ! First level to use

  ! Output
  real(r8), intent(out) :: new_array(:,:,:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = old_array(cols, top_lev:,:)
  end if

end subroutine pack_array_3Dr8

#ifdef GFDL_COMPATIBLE_MICROP
pure subroutine pack_diag_array(old_array, cols,  top_lev, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)          ! List of columns to include
  integer,  intent(in)  :: top_lev          ! First level to use

  ! Output
  real(r8), intent(out) :: new_array(:,:,:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array(:,:,:) = old_array(cols, top_lev:,:)
  end if

end subroutine pack_diag_array
#endif

! Subroutines to unpack arrays for output
!========================================================================
! Rank 1 array of reals, columns only
pure subroutine unpack_array_1Dr8(old_array, cols, fill, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:)   ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  real(r8), intent(in)  :: fill           ! Value with which to fill unused
                                          ! sections of new_array.

  ! Output
  real(r8), intent(out) :: new_array(:)
  
  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = fill

     new_array(cols) = old_array
  end if

end subroutine unpack_array_1Dr8

! Rank 1 array of reals, columns only, "fill" value is an array
pure subroutine unpack_array_1Dr8_arrayfill(old_array, cols, fill, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:)   ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  real(r8), intent(in)  :: fill(:)        ! Value with which to fill unused
                                          ! sections of new_array.

  ! Output
  real(r8), intent(out) :: new_array(:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = fill

     new_array(cols) = old_array
  end if
  
end subroutine unpack_array_1Dr8_arrayfill

! Rank 2 array of reals, columns and levels
pure subroutine unpack_array_2Dr8(old_array, cols, top_lev, fill, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  integer,  intent(in)  :: top_lev        ! First level to use
  real(r8), intent(in)  :: fill           ! Value with which to fill unused
                                          ! sections of new_array.

  ! Output
  real(r8), intent(out) :: new_array(:,:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = fill

     new_array(cols, top_lev:) = old_array
  end if
  
end subroutine unpack_array_2Dr8
#ifdef GFDL_COMPATIBLE_MICROP
pure subroutine unpack_diag_array(old_array, cols, j, top_lev, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  integer,  intent(in)  :: j
  integer,  intent(in)  :: top_lev        ! First level to use
! real(r8), intent(in)  :: fill           ! Value with which to fill unused
                                          ! sections of new_array.

  ! Output
  real(r8), intent(out) :: new_array(:,:,:)

  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
!    new_array = fill

     new_array(cols, top_lev:, :) = old_array(:,:,:)
  end if
  
end subroutine unpack_diag_array
#endif

! Rank 2 array of reals, columns and levels, "fill" value is an array
pure subroutine unpack_array_2Dr8_arrayfill(old_array, cols, top_lev, fill, new_array)
  ! Inputs
  real(r8), intent(in)  :: old_array(:,:) ! Array to be packed
  integer,  intent(in)  :: cols(:)        ! List of columns to include
  integer,  intent(in)  :: top_lev        ! First level to use
  real(r8), intent(in)  :: fill(:,:)      ! Value with which to fill unused
                                          ! sections of new_array.

  ! Output
  real(r8), intent(out) :: new_array(:,:)
  
  ! Attempt to speed up packing if it is unnecessary.
  if (size(new_array) == size(old_array)) then
     new_array = old_array
  else
     new_array = fill

     new_array(cols, top_lev:) = old_array
  end if
  
end subroutine unpack_array_2Dr8_arrayfill

#ifdef GFDL_COMPATIBLE_MICROP
end module micro_mg_mod
#else
end module micro_mg1_5
#endif
